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ABSTRACT 

We study bisymmetric modes of angular wavenumber 2 for flat stellar disks in potentials with smooth 
cores. Stars either all circulate in the same direction or a small fraction may counter-rotate. The 
bisymmetric modes are unstable unless there is a sufficiently large external halo or bulge. We find 
two modes for each disk: a more central fundamental mode and a more extensive and more spiral 
(trailing) secondary mode. The fundamental mode is particularly sensitive to the orbital population 
in the central part of the disk. Depending on that population, it varies from a small compact bar to a 
trailing spiral that is almost as wound as the secondary mode. All modes rotate too rapidly for there 
to be an inner Lindblad resonance. They transfer angular momentum from the central to the outer 
regions of the disk. Most of them release gravitational energy and convert it to kinetic energy, which 
also flows outwards through the disk. 

Subject headings: stellar dynamics, galaxies: kinematics and dynamics, galaxies: spiral, galaxies: 
structure 


1. INTRODUCTION 

The publication by Lin & Shu (1964) of a density wave 
theory of spiral structure generated wide interest, and 
has led to a broad literature. One strand of that litera¬ 
ture consists of stellar dynamic studies of the instabili¬ 
ties of flat stellar systems. Important ingredients of the 
analysis of Lin & Shu are the approximations of tight- 
winding, or near-axisymmetry, of the waves, and of the 
near-circularity of the stellar orbits of the unperturbed 
disk. A key result obtained using these approximations 
is that of Toomre (1964) who showed that the stability of 
stellar disks to axisymmetric waves requires that the pa¬ 
rameter Q = > 1. Here k is the epicyclic 

frequency, an is the radial velocity dispersion, Ed is the 
disk density, and G is the gravitational constant. This 
result has proved to be remarkably robust and widely ap¬ 
plicable, despite the simplifying approximations used in 
its derivation. Yet it was soon recognized that Toomre’s 
criterion does not provide a complete stability criterion 
for stellar disks when iV-body experiments (Miller, Pren- 
dergast & Quirk 1970, Hohl 1971, Ostriker & Peebles 
1973) revealed that disks with Q safely greater than 1 
are prone to large-scale barlike instabilities. 

Kalnajs (1971,1977) developed the matrix method for 
using linear perturbation theory to study the stability of 
stellar systems. His method is quite general and does 
not impose any restrictions on the forms of either the 
instability, or the unperturbed stellar orbits. The matrix 
method is sufficiently complex that it has not been widely 
used, though that is gradually changing. The first com¬ 
prehensive study using it was that of Zang (1976) who 
analyzed the instabilities of scale-free singular isothermal 
disks. Some of the results of his PhD thesis are given in 
papers by his supervisor (Toomre 1977, 1981). Zang’s 
methods were re-used and extended by Evans & Read 
(1998a,b) who analyzed singular scale-free disks with 
other power-law rotation laws. Scale-free disks have the 
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property that the orbits at any energy are simply scaled 
versions of those at any other energy. Though this is 
indeed a simplification, even applications of the matrix 
method to scale-free disks are by no means simple. 

We study stellar disks with smooth and non-singular 
cores. One of the main purposes of our work is to un¬ 
derstand the influence that orbital populations have on 
the responses of the disks to density perturbations, by 
studying the responses of a variety of disks. That vari¬ 
ety allow us to study the effects of the mass and extent of 
the disk, the orbital population, and a central bulge. In 
EH we use the matrix method to derive the eigenvalue 
problem whose outcome is the pattern speed, growth rate 
and mode shape of an unstable galactic disk. We show 
in E21 that a boundary integral, which has been un¬ 
justifiably ignored in some earlier work, may occur with 
unidirectional disks. In >12.31 we extend the linear pertur¬ 
bation theory to second order. Those results are needed 
for studying how modes transfer angular momentum and 
potential and kinetic energy. In El we discuss Poly- 
achenko’s (2004) simplified theory of spiral and bar-like 
modes. In later sections we compare our findings with 
its predictions. 

We introduce the three axisymmetric potentials which 
we use for our models in E They are those of Kuzmin’s 
disk, the isochrone disk, and the cored logarithmic poten¬ 
tial. We use known distribution functions (DFs) for stel¬ 
lar disks with the first two potentials. However, none are 
currently available for cored exponential disks embedded 
in the cored logarithmic potential, and so we construct 
some DFs for them in ^3.51 Lastly in EH we discuss a 
way to model the effects of hot central bulges by modify¬ 
ing DFs by cutting out a part of their orbital population. 

We describe our computational procedures in gl In 
E we give the results of our computations of modes 
of angular wavenumber m = 2 obtained using the ma¬ 
trix method. Some of our models are new, while others 
are the same as those studied in prior work. The prior 
work, some of which was done using carefully applied N— 
body simulations, includes that of Kalnajs (1978), Earn 
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& Sellwood (1995) and Pichon & Cannon (1997) for the 
isochrone potential, and that of Athanassoula & Sell- 
wood (1986), Sellwood & Athanassoula (1986), Hunter 
(1992), and Pichon & Cannon (1997) for Kuzmin’s po¬ 
tential. Vauterin & Dejonghe (1996) studied mod es of 
a cored exponential disk, like that we use in El in a 
potential with a nearly flat rotation curve, which they 
obtained by combining two Kuzmin potentials. Sellwood 
(1989) studied modes of the uncored and mildly singular 
exponential disk in the same cored logarithmic potential 
as we use for our models in El 

In (JHl we discuss and interpret our computational re¬ 
sults. We summarize our results in E 

2. PERTURBATION THEORY 

This section presents the dynamical theory on which 
our work is based. El carries out an Eulerian lin¬ 
ear perturbation analysis of the collisionless Boltzmann 
equation in action-angle variables. We derive a matrix 
eigenvalue problem for finding modes, following Kalnajs 
(1971,1977). We apply this analysis to unidirectional 
disks in which all stars circulate in the same direction 
in El and show that the matrix formulation must then 
include additional boundary integral terms. We extend 
the perturbation theory to second order in a to the 
extent necessary for discussing the transfer of angular 
momentum, and kinetic and potential energy. In El we 
discuss orbits and their responses, and relate the analysis 
of H2.ll to Polyachenko’s (2004) theory. 

2.1. Linear Perturbation Theory 

We study the stability of a collisionless stellar disk 
composed of a distribution of stars moving in orbits 
in a circularly symmetric potential Vo{R). The unper¬ 
turbed system is described by a DF /o(E, L) (Binney & 
Tremaine 1987) where 

E =]^{vr + vI)+Vq{R), L = Rv^, (1) 

are the energy and angular momentum which remain 
constant along an orbit. The density corresponding to 
the unperturbed DF is 

£d = y /odv, (2) 

where dv denote an element of the two-dimensional ve¬ 
locity space. This density is the one that produces the 
potential Vo{R) only in the fully self-consistent case; in 
our work it often provides only a part of that potential. 
The development of the DF for the perturbed state is 
described by the collisionless Boltzmann equation 

df/dt + [f,n]^o, ( 3 ) 

where [, ] denotes a Poisson bracket and H is the Hamil¬ 
tonian 

'E=^{vji + vl)+V{R,(j),t). (4) 

We expand the DF as / = /o + /i + /2 + • • • and the 
Hamiltonian as Ti = Ho + Vi + Vz + ■■ ■, where V = 
Vb + Ui + V 2 + • ■ ■ is the corresponding expansion of the 
potential. Collecting terms of first and second orders 
then gives the equations 

^ + [fi,no] = -[fo,v,], (5) 

^ + [f2:no]=-[fo,V2]-[fuV,]. (6) 


It is necessary that the perturbed densities due to the 
changes from the unperturbed DF are precisely those 
needed to produce the corresponding components of the 
perturbed density, so that 


Uj(x, t) = 

;> 0 . 

( 7 ) 


J |x-x'| 


II 

W 

[ fj{x,v,t)dv, j>0. 

(8) 


Here dx and dv denote elements of the two-dimensional 
position and velocity spaces. Equation o is true when 
j = 0 because we shall use Eq to denote the self- 
consistent density for the potential Vo- Equation m 
replaces the j = 0 case of equation © . 

It is convenient to work with the action variables of 
the unperturbed motion, defined by 

Jr = — (f vndR, = L. (9) 


The first equation here defines Jr as a function of E and 
L. The actions provide an alternative pair of integrals 
of motion, in terms of which we can express the unper¬ 
turbed DF fo{JR, Jtj>) and Hamiltonian Hq[Jr, J^). The 
advantage of using action variables is that their conjugate 
angle variables {Or, 9^) increase uniformly with time 


dOR 

dt 


dHo 

dJR 


^RiJR, J 4 ,), 


dO^ 

dt 


dHo 

dJ^ 


^4>^Jr-! Jtj)')- 


(19 

Perturbations must be periodic in the angles, and this 
allows us to use Fourier series in these angles for them 
(Kalnajs 1971). It is convenient to write these Fourier 
series in the complex form 


/i = ^ MJr,J^)E^^-, (11) 

I —— 00 
00 

^ U(J^, (12) 

I—— 00 


with the understanding that it is their real parts which 
give the physical solution. Following Landau (1946), we 
suppose that the frequency u is complex with real and 
imaginary parts 

uj = iriLlp + is, (13) 


with Lip representing an angular pattern speed and s a 
growth rate. We suppose initially that s > 0 so that 
we have a growing disturbance which was infinitesimally 
small infinitely long ago. The possibility of stationary 
oscillations and real values of uj has to be considered via 
analytical continuation to s = 0 from s > 0. 

The solutions Hill and m are those for a perturba¬ 
tion of a single angular wavenumber m. Perturbations of 
all angular wavenumbers are possible and could be con¬ 
sidered (Kalnajs 1971, Lynden-Bell & Kalnajs 1972) but 
we forgo that generality. Different wavenumbers do not 
interact at the first order when as here, the unperturbed 
state is axisymmetric. In fact we shall be concerned al¬ 
most entirely with the case of m = 2. 

Substituting expansions m and m into equation © 
and matching Fourier coefficients yields the relation 


fl{jR.J<t) 


UjR.J<t) ^ \ 

ILIr + uiLl^ - UJ \ OJr dJcj, J 


(14) 
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The potential Vi, and the density Si which causes it 
and is obtained from integrating fi as in equation 0, 
can also be expanded in position space as 

OO 

(15) 

3=0 

OO 

(16) 

3=0 

where (R) and tr™ (i?) are a complete set of real basis 
functions, and Cj are constant coefficients. We multi¬ 
ply equation ([TT)ll by integrate over 

position space to get where 

OO 

Djkim) = 27ry i3f{R)a'^{R)RAR, (17) 
0 

are the components of a constant matrix D(to). It is 
diagonal if "0™ and u™ form a biorthogonal set. Alter¬ 
natively, we can rewrite Si using its integral and 
then carry out the integration over phase space using ac¬ 
tion and angle variables. This requires that we calculate 
Fourier coefficients of the basis potential functions 


0 

X cos[Wn + m{9^ - (j>)]d6ji, (18) 

OO 

yi=T.^^'^T3^ (19) 

3=0 

(Kalnajs 1977; Tremaine & Weinberg 1984). Using equa¬ 
tion Cl to relate the Fourier coefficients of /i to those 
of Vi yields 


[M)??^,^) — D(77i)]c = 0, (20) 

where the components of the matrix M(m, cu) are defined 
as 


47r^ 


1 — — 00 \ 


dJ, 


<1 dfo 
ddjR 


R 




Iflji + mfl^ — to 


^z^Y:kdJ^- 


(21) 

The integration is over the whole of action space [cf equa¬ 
tion ®], and we suppose that the integrand decays suf¬ 
ficiently rapidly as Jr —>■ oo and ^ ±oo for the in¬ 
tegrals to converge. This system of linear equations 
has a non-trivial solution for the coefficient vector c only 
if 


A4(to, w) = |M(m,a;) — D(to)| = 0. (22) 


This determinantal equation defines a nonlinear eigen¬ 
value problem for w. Details of how to solve it are dis¬ 
cussed in 0 Once w is found, its eigenvector c gives the 
physical shape of the perturbation. 


2.2. Boundary Integrals 

The DF of a unidirectional disk for which all the stars 
rotate in the prograde direction has the form 

MJr, J^) = f^iJR, (23) 


where H is the Heaviside function. The derivative of this 
DF with respect to is 


dfo dfP 


H{J^) + fP{jR,0)6{J^). (24) 


dJ^ dJ^ 

The matrix M(m,a;) then has two components 
M(to, lo) = M^(m, Lo) -f- M®(m, w), 
whose elements are 


(25) 




= E 


l= — oo \ 


47r2 

1 — — OC 


dJfi 


IHr + — oj 

IUr + mfltj, — U3 




(27) 




Because = VIr/2 for radial orbits, the denominator of 
the boundary integral reduces to —uj for I = —m/2 
when m is even. Modes with w = 0 are therefore pre¬ 
cluded when Jq {Jr,0) ^ 0. More generally, this special 
component of the boundary integral can be incorporated 
into an iterative scheme for computing eigenvalues, which 
we describe in 00 

The DF (1^ also drops abruptly to zero at the circu¬ 
lar orbit limit Jr = 0 so that one should also include an 
extra H{Jr) factor in the DF. Differentiation of /o with 
respect to = 0 then gives a S{Jr). However it gives 
no boundary integral, regardless of the value of /o at 
Jr = 0. The reason is that its integrand, for which the 
RT-fo {jR,d) of equation lt?7ll is replaced by 1/^(0, J^), 
vanishes at J/j = 0. That is because the Fourier coeffi¬ 
cients J^) vanish for alH ^ 0 because orbits with 

Jr — 0 are circular, and the remaining non-zero I = 0 
Fourier coefficient is annulled by the factor 1. Hence the 
simpler form of the DF suffices. 

The boundary integral 123 does not arise if the unper¬ 
turbed DF contains no radial orbits so that {Jr, 0) = 
0. This is the case for the unidirectional DFs used by 
Zang (1976) and Evans & Read (1998a,b) because their 
DFs contain positive powers oi L = J^ as factors, and 
hence contain no radial orbits. They work with scalefree 
potentials which are singular as i? —> 0, whereas ours 
are not. Only orbits with low angular momenta pene¬ 
trate near the center of a cored potential. As Gerhard 
(1991) discusses for the analogous problem of spherical 
systems, most DFs for smooth potentials which produce 
finite densities in their cores tend to isotropy there and 
so have radial orbits. His reasoning applies to thin disks 
too. Only orbits with low angular momenta penetrate 
near the center, and the only alternative to radial orbits 
is a singular distribution of non-radial orbits, as when 
all orbits are circular. The need for radial orbits to pro¬ 
vide a non-zero central density disappears if the density 
of the stellar disk drops to zero in the center, as it does 
with the cutouts which we discuss in a 

The omission of the boundary integral terms 123 in the 
cases computed b y Hu nter (1992) invalidates his results. 
As we shall see in 031 the effect is substantial. Pichon & 
Gannon (1997) confirmed those results, but only because 
they repeated Hunter’s error of neglecting the boundary 
integral. As we show in Appendix^ its omission means 
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neglecting the contributions to potential energy and an¬ 
gular momentum which arise from the perturbation of 
radial orbits. 

Boundary integrals are avoided if one uses a La- 
grangian instead of an Eulerian perturbation theory 
(Kalnajs 1977). The two forms of the theory are comple¬ 
mentary and we show in Appendix^that the Lagrangian 
results relevant to this work can be derived simply from 
the Eulerian theory. 


the physical parts of our solutions are given by the real 
parts of our complex solutions. The real parts of /i and 
Vi are given by the sums ^(fi + fi) and ^(Ui + Ui), where 
a bar denotes a complex conjugate. Therefore 

W2,i = ^JJ{h + /i)(Ui + Ui)dJd0, 

= lJJ{fiVi + vJi)d3de, 


2.3. Angular Momentum and Energy 

The perturbations to angular momentum and energy 
are of second order, and so their calculation requires con¬ 
sidering the second order terms from ®. Summing the 
contributions from all elements in phase space gives 

C = JJ J^f dJdQ, (28) 

for the total angular momentum, and 

IC^^JJivl + vDfdJde, (29) 

for the total kinetic energy. To compute the gravitational 
energy, we must distinguish between the part of the 
unperturbed gravitational potential Vq which is due to 
the stars of the disk, and the remainder which is 
provided by some other and external source. The per- 
turbational terms Vj, j > 0, of the potential all arise 
from the perturbed DF, and so all belong to the self- 
gravitational potential. The double contribution of the 
external potential to the gravitational energy is taken 
care of by writing it as 

W = ^ JJ{V + Uo'="‘)/dJd0. (30) 

The leading corrections to these quantities due to the 
perturbations are of second order because all the first 
order terms vanish when integrated over They are 

C2 = JJ Jc^f2dJdQ, (31) 

IC 2 =JJ{no- Uo)/2dJd0 = JC 2,1 + JC 2 , 2 , (32) 

after using the zeroth order part of equation 0 to sub¬ 
stitute for the velocities, and, because kg®’'* = Vo — Vjf, 

W2 = ^ + Vih - Vo^f 2 + 2Uo/2)dJd0, 

= l JJ fiVidJde + JJ Vo f 2 dMQ, 

= }V2A+m,2. (33) 

The component VV 2,2 = —1^2,2 because their defining in¬ 
tegrals match. The step to the second line of equation 
(|S3 uses the fact that the first and third terms on the 
first line cancel. This is seen by transforming the integra¬ 
tions to (x, v) space, and then using equations 0, 0 
for j = 2, and Poisson integrals like o to relate poten¬ 
tials to densities, to express them as identical integrals 
of the product of the densities Eq and S 2 . 

The simplest integral to calculate is that for the first 
component of W 2 ,i in equation 03- As we noted earlier, 


CXD 

= ^ (34) 

I — — OC1 


where the components W 2 i are defined by 
WL=Tr'^ f d-T (^ 




X 


[^Ufl + m(U0-Up)]|Up 




(35) 


We have used here, and will again, the fact that the 
only products which do not vanish on integration over 
9(f, are those which pair a conjugate with a non-conjugate 
quantity. 

We show in appendix IXI that 


00 00 

C2{t)=e^^* 4, /C2.i(t)=e2*‘ ^ (36) 

I —— 00 I —— 00 


where the components of these sums are defined by the 
integrals 


L 2 = —TOTT^ j dJ 

_ ' ' 

\lflR + mfl^ 


dfo dfo\ 


iC' ^ = -TT^ / dJ ( Z 


dfo 

dJR 



{IUr + mil^)\Vi\'^ 
\lflR + mfl^ — tuP 


(37) 


(38) 


To each of the area integrals m . (|23 and (123 must be 
added the boundary integrals given by the delta function 
term of equation (EH) for the prograde DF (E3. The 
integrals can be combined to give the simple relation 

Kh + Wl, = npLi, (39) 


between the separate components. Although each of 
these components can be found directly from the first 
order solution, computing VV 2.2 = —^ 2,2 requires more, 
as we show in Appendix 0 
The second order corrections £2 and >V 2 ,i to the an¬ 
gular momentum and gravitational energy have simple 
representations in terms of the real and imaginary parts 
of the matrix M = Mr + iMj. Combining equations 
E3 and (1^ with expansion 113, we get the quadratic 
forms 

£2(t) = -^e^'**c^Mic, (40) 

4s 

W2,i(<) = ^e^®*c^MRC, (41) 
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Fig. 1. — T he orbital frequency space of the cored logarithmic 
potential isni with frequencies in units of vq/Rq. The lines are 
curves of constant E in equal steps in e~^. They show that 
depends only weakly on the angular momentum L, which varies 
from zero at their lower limit on the boundary = 0.5^2/^ which 
corresponds to radial orbits, to Lc{E) on the curved upper limit 
which corresponds to circular orbits. The largest value achieved by 
$1^ = — 0.5r2/^ is 0.106. Similarly narrow lens-shaped plots are 

obtained for other cored potentials. When normalized to the same 
range of frequencies as in FigE the largest value of is 0.130 for 
Kuzmin’s disk and 0.119 for the isochrone. depends only on 
E for the isochrone, so that curves of constant E are then exactly 
straight. 


where the superscript T denotes transposition of the col¬ 
umn vector c to generate a row vector. These expressions 
are real because the matrix M is symmetric. Moreover 
multiplying equation (EOl) by and separating real and 
imaginary parts (D is also symmetric) gives 

c^MrC = c^Dc, c^MrC = 0. (42) 

The second relation shows that C 2 {t) = 0, which it 
must be because the disk is not subject to any exter¬ 
nal torques, and hence its angular momentum is con¬ 
served. There is no such restriction on sizes of the sep¬ 
arate Fourier components represented by the terms for 
different I in the sum PHl . other than that they must 
sum to zero. Similarly the sizes of the different compo¬ 
nents of the potential and kinetic energy vary, because 
only the total energy £ is constrained to be zero, with 
£ 2 {t) = K. 2 {t) + 14^2(0 = £lpC 2 (t) = 0. Hence the two 
components of the kinetic and potential energy are re¬ 
lated in the same way: 

/C 24 (t) = ->V2,l(t), /C2.2(t) = -W2.2(t). (43) 

The fact that the second order components grow twice 
as fast as those of first order is not paradoxical. It re¬ 
flects the fact that our analysis can describe only the 
early stages of the growth of an instability. If we intro¬ 
duce a small ordering parameter e into our expansion 
f = fo + ^/i + e ^/2 + • ■ ■ of the DF to measure the size 
of the perturbation relative to that of the unperturbed 
state, then we see that the linearization breaks down, 
and our analysis is unreliable, after a time t such that 
ee'’* = 0(1). The second order components are then of 
magnitude which is also 0(1). Our perturbation 

theory has then ceased to be valid. It is useful only when 
our expansion is well-ordered, that is for times for which 
ee'** is small. 

2.4. Orbits and Polyachenko’s Unified Theory 

Every orbit in a disk with an axisymmetric potential 
is a rosette (Binney & Tremaine 1987). One extreme 
case is that of circular orbits with the maximum angu¬ 
lar momentum Lc{E) for the energy E. They have no 


radial motion and their radial action Jr = 0. The other 
extreme is that of radial orbits with zero angular momen¬ 
tum J^. For them the polar angle f remains constant, 
except for discontinuous changes by tt whenever the orbit 
passes through the center. Intermediate orbits oscillate 
in R between maximum and minimum orbital radii i?max 
and Figure n shows an orbital frequency space for 

prograde orbits. Its shape is characteristic of those of 
other cored potentials. The upper curved boundary is 
formed by circular orbits for which = VfR)/R, where 
Vc{R) is the circular velocity at radius i?, and ULr = k{R) 
is the epicyclic frequency. The straight lower boundary is 
formed by radial orbits for which = £Ir/2. The slim¬ 
ness of the region in between these boundaries, where 
all the intermediate orbits lie, shows how limited is the 
range of — £Ir/2 for all orbits. 

The frequency = 12^ — £Ir/2 and its narrow 
range has long been recognized as dynamically impor¬ 
tant (Lindblad 1959). Lynden-Bell (1979) pointed out 
that the response of an orbit to a weak bar-like poten¬ 
tial, which rotates with a pattern speed Op for which 
jHp — Hil is small, is to align the orbit with the bar if 
Hi decreases as decreases when Jf = Jr-\- J^/2, the 
action corresponding to the frequency Oi, is fixed. This 
is the case when 




dOj 

dL ^ * dE 


> 0 . 


(44) 


Lynden-Bell labeled such orbits as abnormal, as opposed 
to normal or donkey orbits which respond contrarily to 
a bar-like force. It is straightforward to classify orbits in 
any potential according to the criterion ijHJ; abnormal 
orbits are those which lie below the dashed lines in Fig¬ 
ure ^ All sufficiently radial orbits are abnormal, but so 
too are all orbits which remain in the central regions of 
cored potentials. The dividing curve for the cored log¬ 
arithmic potential in Figure 31 tends to the asymptote 
L/Le = 0.723 as E ^ 00 (See fi21l. That asymptote, 
which is barely evident due to the compressed scale at 
high energies, is the dividing curve between normal and 
abnormal orbits for the singular logarithmic potential at 
all energies. 

Polyachenko (2004) has outlined a unihed theory of spi¬ 
ral and bar-like modes of disk galaxies based on Lynden- 
Bell’s (1979) analysis. The analytical basis of his theory 
is a simplihed version of the analysis given in 0 He 
uses the smallness of [Dp — Hij to deduce that, to leading 
order, the perturbed DF is given by a single I = —m/2 
component of the Fourier series dJ so that, in our no¬ 
tation, 

Polyachenko’s equation (11) relates the perturbed DF 
to its potential when averaged over the ‘fast’ angle Or. 
That averaged potential is given in our notation by the 
single I = —m/2 component of the Fourier series HI 211 . 
Polyachenko’s equation (11) is therefore equivalent to the 
relation 

(mn. - u,)U,, = “ (|| - («) 

He has a different method for solving his equation (11) 
which has the advantage that he is able to derive a lin¬ 
ear eigenvalue problem for w, albeit one which must be 
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solved for an unknown function of the two action vari¬ 
ables. Equation igni) shows that his method is equivalent 
to the Kalnajs matrix formulation when that formulation 
is simplified by replacing the sum over I in the coefficient 
equation m by the single I = —to /2 term, or just the 
/ = — 1 term in the important to = 2 case. Hence the 
matrix method can be applied to check the accuracy and 
validity of Polyachenko’s approximations by comparing 
its results with those obtained when only the I = —1 
terms of the matrix M(2, oj) are used. Such comparisons 
are made throughout and discussed in io 

3. MODELS FOR GALACTIC DISKS 

In >13.II >13.21 and a we give the three poten¬ 
tials which we use for stability calculations; Kuzmin’s 
disk, the isochrone, and the cored logarithmic potential 
(Kuzmin 1956, Binney & Tremaine 1987). The rota¬ 
tion curves for both Kuzmin’s disk and the isochrone 
rises at small radii, peaks at a characteristic radius, and 
then falls like a Keplerian potential at large distances. 
That for the cored logarithmic potential gives a flat rota¬ 
tion curve at large distances. All potentials have smooth 
cores. In El we discuss constraints on the sizes of ex¬ 
ponential stellar disks which lie in cored logarithmic po¬ 
tentials. In ^3.51 we construct DFs for such stellar disks. 
In >13.61 we discuss cutout functions which can be applied 
to DFs to model hot central bulges. 


3.3. The Cored Logarithmic Potential 

The cored logarithmic potential, which has been widely 
used in galactic studies because of its flat rotation curve, 
is 


Vo{R)=vl\n^l + R^/Rl, 


(50) 

where vg is the flat rotation velocity, and Rc is again 
the core radius. The disk density needed to produce the 
potential is 


So(ii) = 


-.E 


R^ 






(51) 


2TTGy^W+I^ 

where E is the complete elliptic integral of the second 
kind. This self-consistent density may be derived by 
Toomre’s (1963) Bessel function method, using in turn 
formulas (6.566.2), (6.576.3), (9.131.1), and (8.114.1) of 
Gradshteyn & Ryzhik (1980, hereafter GR). It can also 
be obtained by taking the /3 ^ 0 limit in equations (4.7) 
of Qian (1992) for to = 0 and 71 = 0.5. The limit is 
straightforward for the density, but it is necessary first 
to subtract the constant value of the potential a,t R — 0 
before taking the f3 ^ 0 limit for the potential. 

The potential-density pair and inj is similar to 
that of the Rybicki disk (Zang 1976; Evans & Gollett 
1993) for which 

Uo(i?)=i^o'ln(^l + ^l + RV^c) > (52) 


3.1. Kuzmin’s Disk 

Kuzmin’s disk is known also as the Kuzmin-Toomre 
disk because it is Model 1 of a family given by Toomre 
(1963). Its potential and self-consistent density are 


Vg{R) 


-GM 


So(i?) 


MRc 

27r(R2 + )3/2 


, (47) 


where Rq is the core radius of the potential. Kuzmin’s 
disk has been widely used for stability studies because 
of its simplicity. Sellwood & Athanassoula (1986) and 
Athanassoula & Sellwood (1986) used it for A'-body sim¬ 
ulations, while Hunter (1992), Pichon & Gannon (1997) 
and Polyachenko (2004) have used it for analyses based 
on the collisionless Boltzmann equation. To have results 
which can be compared with those of earlier work, we 
use the DFs given by Miyamoto (1971) and those used 
in Athanassoula & Sellwood’s simulations. 


3.2. The Isochrone Disk 


The potential and self-consistent density of the 
isochrone disk are 


Vo{R) 

So(i?) 


-GM 

Rc + y/WT^' 

MRc\ ^R+^/WT^ 

2^R3^ 


(48) 


R 

x/i ?2 + 


(.49) 


Again Rc is the core radius of the potential. We use the 
DFs given in Kalnajs (1976b), modified so as to allow for 
a population of retrograde stars in the manner specified 
in Earn & Sellwood (1995). The stability of these models 
was investigated by Kalnajs (1978), Earn & Sellwood 
(1995) and Pichon & Gannon (1997). 


So(R) = - 


(53) 

2t:G^R^ + R^c 
Rybicki’s disk is obtained by subtracting Toomre’s 
(1963) Model 0 from a singular Mestel (1963) disk. The 
density EU for the cored logarithmic potential is larger 
than that of the Rybicki disk in the center, which is why 
it has somewhat larger circular velocities. Both disks re¬ 
semble the singular Mestel disk, whose stability has been 
studied by Zang (1976) and Evans & Read (1998a,b), at 
large distances. Both tend to that singular limit as the 
core radius Rc 0 . 


3.4. Exponential Disks in the Cored Logarithmic 
Potential 


It is widely accepted that the densities of the stellar 
disks of spiral galaxies decay exponentially with distance 
(Freeman 1970). We study cored disks with densities of 
the form 


T,d{R) = T,s exp 


VR^ + 77^ 

Rd 


(54) 


where Ru is the length-scale of the exponential decay, 
and is a density scale. We include the core radius 
Rc of the potential in the density to avoid the 
logarithmic singularity of the uncored exponential disk 
Ed = Eg exp(—R/R d) (Freeman 1970). Vauterin & De- 
jonghe (1996) did likewise. 

The self-gravitational potential due to the cored disk 
density can be derived by the method of Evans & 
de Zeeuw (1992) as 


V^{R) = -2ttGEsRd 


^Ji{x)da 




A 2 ^a ;2 -fA 2 


(55) 
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Here A = Rc/Rd and Ji is the Bessel function of the 
first kind. Beware that this integral must be evaluated 
with great care because it is oscillatory and its amplitude 
decays only as as x —> oo. V(f(R) reduces to the 

potential of the exponential disk in the limit of Rc —*■ 0 
when A ^ 0. 

There are limits on the scales of a cored exponential 
disk lying in the cored logarithmic potential CT . 
The difference [Vb(.R) — between the total poten¬ 

tial and that of the disk must be provided by some halo 
or bulge component of the galaxy, whose density must 
be everywhere positive. If that halo/bulge is spherical, 
then its density is given by Poisson’s equation as 


PH(r) 


1 1 d 

47 rG dr 




v^(R) = R 


dVo(R) 

dR 


dVo^jR) ' 

dR 


> 0 , 


(56) 

(57) 


(Zang 1976; Vauterin & Dejonghe 1996) where r = 
VR^ + Here v^{R) gives the amount by which the 

square of the circular velocity for the galaxy exceeds that 
of the disk alone. The physically necessary condition 
that pH(r) > 0 requires that v‘^{R) > 0, but is more 
restrictive. The region below the solid curve in Figure 
m gives the range of the dimensionless combinations of 
parameters CSsRd/vq and Rq/Rd which it allows. The 
boundary value GT,sRb/vq = 0.304 of the solid curve 
at A = 0 applies to the limit of the uncored exponential 
disk in a singular logarithmic potential. An alternative 
statement of the condition pH(r) > 0 in this limit is 
^o(^d/G-^d)^^^ > -723, where Md is the mass of the 
uncored exponential disk. 

Sellwood (1989), and unpublished work by Toomre de¬ 
scribed there, studied modes of the uncored exponen¬ 
tial disk in the cored logarithmic potential iP|l. The 
critical value of voiRD/GMi))^/'^ needed to ensure that 
Ph > 0 for this case is not greatly changed when the 
logarithmic potential is cored. It is reduced only to 
0.691 when A = Rc/Rd = 0.5. For A = 0.2 and 
vo{Rd/GMdY/^ = 0.6 as in the model displayed in Sell- 
wood’s Figures 1 and 2, not only does pn become neg¬ 
ative, but x^(i?) < 0 for 1.4 < R/Rd < 3.0. The fixed 
halo therefore exerts an outward force in this range. (The 
singular nature of the potential of the uncored exponen¬ 
tial disk also causes pn < 0 , but only for R/Ru ^ 1 
which is well below the softening length used in compu¬ 
tations.) 

Giovanelli & Haynes (2002) analyzed a large number 
of rotation curves and have found that the ratio of the 
scale length of the steep inner rise of the rotation curve 
to the scale length of the exponential varies from 0.63 
for the most luminous galaxies, to 1.28 for their least 
luminous. They fit a parametric model in which the ro¬ 
tation curve decays exponentially towards its form in the 
outer regions, and hence their ratios provide only the ap¬ 
proximate guidance that it is to reasonable to use values 
of order unity for our ratio Rc/Rd- Figure El plots to¬ 
tal circular velocity, and the parts provided by the disk 
and the halo/bulge for three different relative sizes of 
the exponential disk which more than cover their range. 
The value of GUsRd/vq is close to 90% of the maxi¬ 
mum allowable for that Rc/Rd in each case. The least 
extensive disk makes the largest contribution to the ro- 



Fig. 2. — The boundary of physical exponential disks (solid line) 
in the parameter space. Models below the solid line have spherical 
halos with pH > 0 at all radii. 



Fig. 3. — The rotational velocities due to the disk (upper panel) 
and halo/bulge (lower panel) components of the exponential disk, 
in units in which Rq = uq = 1. In each case the dashed, dot- 
dashed, and dotted lines correspond to (Ry)^GT>sRd) = (0.5, 0.6), 
(1,0.42) and (2,0.32), respectively, while the the solid line shows 
the circular velocity of the cored logarithmic potential. 


tation curve at the center, but its relative contribution 
then declines rapidly. The disk’s contribution tracks the 
rotation curve considerably further in the intermediate 
case with Rd = Rq^ and requires only a relatively small 
contribution from a central halo/bulge to make up the 
deficit. As i?D increases, the disk becomes less maximal 
and an increasingly large spherical central halo/bulge is 
needed. 


3.5. Distribution Functions for the Exponential Disk 
We construct DFs by using the identity 




7?c ^0 


(58) 


to partition the density in the form 
SD(i?) = E,e-^^%-2^^ (^1 + ^^ , A=|^. (59) 
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Fig. 4.— (a) The DF /q of the exponential disk for N = 6 and Rq = i?D = 1, and (b) that DF after applying the cutout insj with 
Lq = 0.1. The abscissa u = y/l — e~^ is that used in numerical work, as defined in equation EJ. It ranges from u = 0 at the center to 
u = 1 at infinite distances. Contours are labeled in units of Ss, and show slopes which increase with increasing N. The dashed line shows 
the curve d^ijdL + QidQi/dE = 0. Only orbits which lie above this curve are normal in Lynden-Bell’s (1979) classification. 




Fig. 5. — (a) Circular velocity (solid line) and mean rotation velocities of the exponential disk. Dashed, dot-dashed, dotted and 
long-dashed lines respectively show the mean rotation velocity of exponential disks with = 2, 4, 6 and 8. All curves are for 

Rq = Rd = uq = 1. (b) Toomre’s Q for the exponential disk for Rq = 1 with N = 2 (thin lines) and N = 6 (thick lines). Solid, dashed 
and dot-dashed lines here correspond to (i?D?GSsi?D) = (0.5, 0.6), (1,0.42) and (2,0.32), respectively as in Figure!^ 


Here N is an integer parameter which allows us to 
generate a family of models. Binomial expansion of 
(1 + B?/Rq)^ gives Ed as a sum of powers of R? mul¬ 
tiplied by functions of the potential <&. We then use 
Sawamura’s (1988) method, following Evans & Collett 
(1993) §3.1, to find the DF as the series 

fP{E, L) = S, Aj g„(E), (60) 

where the functions gn{E) are given by 

“ 2"0Fr(n-f 1/2) 


d"+i 




'”)■ ( 61 ) 


dE^+i 

Figure^ plots an V = 6 case of equation CT . 

FigurelSt plots the mean rotation velocity {v^) for four 
different values of TV, and shows that the disks become in¬ 
creasingly cool and centrifugally supported with increas¬ 
ing TV. The mean rotation velocity exceeds the circular 
velocity near the center, due to the n = 0 isotropic term 
in the DF (16011 . Its density, which is 

1 -I- I exp 




i?2 


i?D 


(62) 


is confined to the central region. It is compressed as N 
increases, and disappears as TV —> oo along with all radial 
orbits. The other interesting feature of Figure Eh is that 
(v^) declines away from vq at large distances. Evans & 
Collett’s (1993) models for a simple uncored exponential 
disk in a singular logarithmic potential have the same 
property - see their Figure 2a - and for the same reason. 
Although (u^) —!■ Vo as N —!■ oo, their equation (3.17) 
also shows that (v^j decreases for increasing R for fixed 
TV, as in their Figure 2a and our Figure Eh- 

Figure Eh plots the Toomre stability parameter Q for 
two different TV values and for the same three sizes of the 
exponential disk as in FigureEl The cooler TV = 6 models 
have regions of varying extent in which Q is close to its 
marginal value of 1. The growth of Q with increasing 
R is due primarily to decreasing Ed because the radial 
velocity dispersion cr/j = falls off only mildly. 

3.6. Simulating a Hot Bulge 

Most spiral galaxies of S type have a three-dimensional 
central bulge which stands out from the disk. Kor- 
mendy (1977) proposed modeling the disk with an inner- 
truncated exponential so that the disk has a hole in its 
center, and others have followed. We construct disks 
with central holes by removing stars which penetrate to 



























UNSTABLE BAR AND SPIRAL MODES OF DISK GALAXIES 


9 


the center. We use a cutout function Hcut{L) to change 
an unperturbed DF fo{E,L) to iLcut(L)/o(A, A). In so 
doing, we make the tacit assumption that the part of the 
gravitational potential which the cutout stars had previ¬ 
ously provided, is instead provided by a bulge which is 
so hot that it does not respond at all to disturbances in 
the disk. 

Cutouts were introduced first in Zang’s (1976) pioneer¬ 
ing study of DFs in a singular scalefree logarithmic po¬ 
tential. His purpose was to eliminate stars with short 
dynamical time scales and high dynamical frequencies. 
Dynamical frequencies in cored potentials are bounded 
(cf Figure ^1 and so avoid that problem. Zang’s cutout 
also removes stars on nearly radial orbits, regardless of 
their dynamical frequencies. Evans & Read (1998a) went 
further and also used an outer cutout to remove stars 
that spend so much time far out in the disk that they 
do not respond to its changing potential. Their disks are 
infinitely massive whereas ours are not; ours have suffi¬ 
ciently few such stars that they do not pose a problem. 

We apply a cutout function 

i7ent(A) = l-e-(^/^«)', (63) 

where Lq is some angular momentum scale. Like Zang’s, 
it removes stars for which L is significantly less than Lq, 
but has no effect on stars with L ^ Lq. It also removes 
stars on radial and near-radial orbits, as well as stars of 
low energy because they too have low angular momenta. 
The result of the cutout lIMll is to give an active surface 
density 

Sact = y’i7cut(A)/o(F,L)dv, (64) 

which tends to zero at the center, and hence models a 
central hole. Figure @1) shows the effect of an Lq = 0.1 
cutout on the DF of Figure No boundary integrals 
m arise for such cutout unidirectional disks because 
their unperturbed DF vanishes at A = 0. The active 
surface density corresponding to the cutout DF of Figure 
H) is shown in Figure ITlfc. It shows that our truncation 
is far less sharp than that given by Kormendy’s (1977) 
inner-truncated exponential formula. 

4. COMPUTATIONAL PROCEDURES 

This section gives details on the numerical methods 
and algorithms we used. Those uninterested in these 
topics should skip to where we present the results of 
our computations. 


4.1. Choice of Basis Functions 

We use the set of basis functions introduced by 
Glutton-Brock (1972), and simplified by Aoki & lye 
(1978), to relate densities and potentials. They are 


=_p™ (C) 

{l + R^/b^C 


(2m -I- 2j + 1) 


pn 


27rG&(l + i?2/&2)3/2"m+j 


(6, 


(65) 

( 66 ) 


where ^ = (i?^ — b‘^)/{R^ + 6^) and the Pm+jiO are as¬ 
sociated Legendre functions. The variable ^ ranges from 
— 1 at i? = 0 to 1 as i? —> oo. The functions (in3 and 
(1^ form a complete biorthogonal set over the range 
0 < i? < oo, and so are useful for disks of infinite extent. 


The parameter 6 is a length scale which can be chosen 
so as to match the oscillatory behavior of the Legendre 
functions with the regions in which modes vary most. 
We have found that repeating calculations with different 
b values provides a helpful check on their accuracy. 

The Glutton-Brock functions with the choice b = Rq 
are a simple choice to use with Kuzmin’s disk. The 
self-consistent potential-density pair EH) is then given 
by (Vb,So) = [GM/Rc){iI)q,(Tq). However we have also 
found Glutton-Brock functions to work well for determin¬ 
ing the modes of the exponential disk ini), even though 
its density decays more rapidly than the R~^ decay of 
the cr° functions. This seems to be because modes de¬ 
cay rapidly in the outer regions and are little influenced 
by matter there. For instance, the frequencies which Pi- 
chon & Cannon (1997) computed for a Kuzmin’s disk 
truncated at i? = 5i?c differ little from those computed 
by Hunter (1992) who did not truncate them. Unfor¬ 
tunately, the Glutton-Brock functions are not suitable 
for computing reliable values of the gravitational energy 
component >V 2,2 using the sum 1A18II . The reason is that 
the combination of the growth of the logarithmic poten¬ 
tial Vq at large R with the R~^ decay of the cr° causes 
the integrals in that sum to grow slowly with increasing 
j. The convergence of the sum is poor, and we do not 
have accurate values for >V 2 , 2 - No such problem arises 
with the more rapidly decaying potentials of the Kuzmin 
and isochrone disks. 

Kalnajs (1976a) gave a set of biorthogonal Abel-Jacobi 
functions which are suitable for disks of finite radius, 
and were used by Earn & Sellwood (1995) and Pichon & 
Cannon (1997). No recursive relations have been given 
for them, and so they must be computed using explicit 
formulae. That requires high-precision arithmetic be¬ 
cause of cancellations between large coefficients of oppo¬ 
site sign. Glutton-Brock functions are computed accu¬ 
rately using simple recursive formulae, and do not need 
high-precision arithmetic. 

4.2. Action Space Integrations 

Since DFs are usually given in terms of the energy E 
and the angular momentum A, it is more convenient to 
use E and A as the integration variables for integrations 
over action space. We do this and replace the Jacobian 
dJjidJ^ with dEdL/Llji. To evaluate the surface inte¬ 
gral in the {E, A)-space, we adopt the trapezoidal rule 
in the Al-direction and an extended open scheme [equa¬ 
tion (4.1.18) of Press et al. (1992)] in the A-direction. 
With this latter choice, we avoid the circular and radial 
orbit boundaries where the calculation of orbital frequen¬ 
cies and their derivatives with respect to the actions has 
some computational difficulties. The Fourier coefficients 
dtY), and hence the integrands, all vanish at Amin and 
Amax because the potential functions if^{R) vanish at 
A = 0 (for TO > 0) and as A ^ oo. For the exponential 
disk for which 0 < A < oo, we transform from E to u 
where 

E = (67) 

and map 0 < A < oo onto the finite range 0 < m < 1. We 
generate a uniform no x me mesh in either (A, A/Ac)- 
space or (it, A/Ac)-space, and evaluate integrands at grid 
points. 
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Accurate evaluation of the boundary integrals M® re¬ 
quires an especially careful treatment of the central re¬ 
gion, especially for centrally concentrated modes. More 
basis functions are needed there, up to 15 in some cases, 
and a fine grid to allow for the oscillations of Fourier co¬ 
efficients near the center where Jq is largest. We use a 
uniform grid in u = Vl + E for boundary integrals for 
Kuzmin’s disk; the transformation still works well 
for the exponential disk. 

Formulas for the orbital frequencies and de¬ 
fined in m are obtained by differentiating equation 0 . 
Computationally convenient forms of the integrals which 
occur in them are given by Evans & Read (1998a). An 
alternative and compact form is obtained by expressing 
E and L in terms of the maximum and minimum orbital 
radii, i?niax and i?niin- Then if we define x = and 
U{x) = 2 B?Vq{R), we get 


2R (F Vo) L — (^max x){x ^min)f^[^min7 ^max] ■ 

( 68 ) 

Here f7[a:niin) Xmax] denotes a second order divided dif¬ 
ference, as in de Zeeuw & Hunter (1990). Using this 
result, and then the substitution x = x^iacos^ ip + 
a^max sin^ (fi, the orbital frequencies are given by 

TT 

J^) = -r, 
to 

where 

7r/2 


a^^ ^U[Xniin5 a;? a^max] 



n^iJR,J^) = ^, (69) 

to 

dip . 


These integrals are free of singularities, and they can be 
differentiated analytically with respect to E and L to 
obtain the derivatives of the frequencies needed in El 


4.3. Fourier Coefficients 

We use the trapezoidal rule to evaluate the Fourier co¬ 
efficients CHI because it is a highly accurate method for 
integrating periodic functions if a fine enough gridding 
is used (e.g., Davis & Rabinowitz 1984). Closed form 
expressions for R(t) and (j)(t) are known for orbits in the 
planar isochrone (Boccaletti & Pucacco 1996). Other¬ 
wise orbits must be obtained by numerical integration. 
We integrate the equations of motion for a half period, 
starting from the initial conditions t = Or = 0^ = vr = 
(f) = 0, R = i?min, and record the values at equal steps in 
Or for use with the trapezoidal rule. 


4.4. Evaluating Matrix Elements 

The denominator term {I^Ir mCl^ — to) in the inte¬ 
grals EO for the matrix elements Mjk may vanish if uj 
becomes real during the course of the iterative search. 
When it does, the integral becomes singular. We use 
Hunter’s (2002) method to handle such resonances. It is 
designed to ensure that matrix elements are calculated 
according to Landau’s (1946) rule, and also provides an 
efficients way for the repeated evaluations of the matrix 
elements needed in an iterative search for eigenvalues. It 
requires the evaluation of the integrals 


an(j, fc, 0 = 


4(2n -I- l)7r^ 

(^max ^min) 


dfa , dffi 




m 


dJd> 


X ^J^id>liP4w{l)]dJRdJ^, 


(71) 


where r] = I^Ir + and w{l) = 2{r] - r7min)/(??max - 
?7min) — 1. The extreme values of ij depend on l as is seen 
from Figure n Matrix elements are then formed as sums 

oo oo ^ 

Mjk{m,uj)= ^ <-2 ^ a„ (j, fc,/) (5„ [A(0] > . 

1— — 00 K n=0 ) 

(72) 

where X{1) = 2{uj - ??min)/(?ymax - Vmin) “ 1- Beware 
that the sign of the matrix M used here is the opposite 
of that in Hunter (2002). Both and Qn denote the 
usual Legendre functions. The multivalued Qn functions 
must be evaluated for real lu by taking the limit s = 
Im(a;) ^ 0 through positive values. The coefficients a 
are computed and stored for the whole ranges of n, j, k 
and 1. Weinberg (1994) gives an alternative method of 
handling resonances. 


4.5. Searching for Eigenvalues 

We search for eigenvalues oj using the modified Newton 
method of Stoer & Bulirsch [1993, eq. (5.4.1.7)]. We cal¬ 
culate A4 (to, uj) using an LU decomposition as in Press et 
al. (1992). We compute the derivative dM/duj, needed 
by Newton’s method, numerically using central differ¬ 
ences. Whether Newton’s method converges or not, and 
the rapidity with which it converges if it does, depends on 
the choice wq of a starting approximation. We are most 
interested in finding the modes with the largest growth 
rates. The set of all initial guesses of ujq which lead to 
an eigenvalue is the basin of attraction of its mode. Our 
computations show that the cuo-space is dominated by 
the basin of attraction of the fastest growing mode, but 
that there are also regions which lead to other modes. 
We have found two modes for most of our models. 

We have found the following search procedure to be 
useful when there are large boundary integral terms. As 
we noted in E2 the ^ = — 1 component of the boundary 
integral M® of (EH) is 


47r^ 

LJ 


0 




jk: 


(73) 


where Ajk are the components of a positive definite ma¬ 
trix A which is independent of uj. We define the matrices 
B and E as 

B = M+-A, E = A^i-(B-D). (74) 

UJ 

The eigenvalue problem can be recast as 

(E-il)c = 0, (75) 


SO that l/w is an eigenvalue of w-dependent matrix E. 
We find that it is simple to find roots of the reduced 
equation |E| = 0 by Newton’s method because they are 
insensitive to the initial choice of uj. We use a continu¬ 
ation scheme to find eigenvalues of the recast eigenvalue 
problem by introducing a parameter p and defining 
a sequence of problems 


E(m, uj) 



UJ 


(76) 


We proceed by continuation in p from the simple p = 0 
case to the p = 1 case of equation G3, using our solu¬ 
tions at one stage as the initial estimate for the solution 












UNSTABLE BAR AND SPIRAL MODES OF DISK GALAXIES 


11 


at the next stage. We also apply continuation methods 
when we use eigenvalue estimates obtained from matrix 
truncations of one size as initial estimates for larger ma¬ 
trix truncations. 

Getting a converged solution needs some controls on 
our series truncations and integration grids. Using a fine 
grid at the energies corresponding to distant stars, re¬ 
quires too many Fourier coefficients to expand the basis 
functions. In turn, a large Imax (l^minl) increases the cost 
of computations while such an attempt does not capture 
further physics of the problem. In fact, unstable waves 
are confined to the central regions of stellar disks and 
they lose much of their power as they reach the corota¬ 
tion and outer Lindblad resonances. So, large-amplitude 
orbits, which spend much time in outer regions, do not 
respond to density perturbations. This is the reason why 
outer-truncated disks have been used successfully by oth¬ 
ers (Athanassoula & Sellwood 1986, Pichon & Gannon 
1997). We allow disks to be infinite. We start the com¬ 
putation of w with 40 X 40 grids and = 2. We si¬ 
multaneously increase jmax and grid size until the relative 
accuracy in computing oj becomes better than 5 x 10“^. 
We find in most cases that this accuracy is achieved when 
nc « 75 and jmax = 10. An independent check of ac¬ 
curacy is to see whether the condition e~‘^^^UpC 2 = 0 
is satisfied after finding lu and c. In all of our unstable 
models, the normalized value of this quantity remains 
below 5 X 10“^. 

After an eigenvalue uj has been found, we need its 
eigenvector c. A theorem of linear algebra states that, 
when io is an eigenvalue, then any column of the adjoint 
matrix adj [M(m, uj) — D(to)] is an eigenvector associated 
with UJ (Brogan 1990). The computation of the adjoint 
matrix requires the computation of the minor determi¬ 
nants of M(m,a;) — D(m). It can be done rapidly and 
accurately because we work with matrices of relatively 
small size. 

5. M=2 MODES OF DISKS 

The following three subsections present our numerical 
results for the different DFs we have studied. Properties 
of TO = 2 modes are listed in the tables, and selected ones 
are displayed. The final column in each table gives the 
pattern speed Dp computed by Polyachenko’s simplified 
theory described in in No growth rate s is listed for 
this case because it is always 0. 

We have found two modes for most of the models we 
have investigated. Glassifying them is less straightfor¬ 
ward than it is for simpler physical systems, for which 
the fundamental mode has the lowest frequency and 
simplest structure, and modes of successively higher or¬ 
der have higher frequencies and more complex structure. 
The most important instability is that with the highest 
growth rate. Often it also has the simplest structure. 
However, we find that small changes in the orbital popu¬ 
lation have a much greater effect on relative growth rates 
than they do on radial structure, with the result that 
the mode with the simplest structure is not always the 
fastest growing. For that reason we have chosen struc¬ 
ture, rather than growth rate, as our criterion for de¬ 
termining which modes are fundamental. Fundamental 
modes are labeled 1 and secondary modes 2 in our tables. 

The displays include a contour plot of the perturbed 


density 

El = P(R) cos [2(1) -Ujt + -d{R)]. (77) 

It is obtained from the real part of equation CH) after 
writing its i?-dependent part in the form P(i?)e*'^^^^ for 
some real functions P{R) and )}{R). Eigenvectors are 
arbitrary to within a complex constant multiple, with 
the result that the phase i?(i?) is arbitrary to within an 
additive constant, so that modes are oriented arbitrarily. 
As usual, we draw only the contours for positive levels of 
the perturbed density G3; those for negative levels have 
exactly the same pattern, rotated by 90° and occupy the 
blank sectors. The levels of the contours are in steps 
of 10% of the maximum of Si from 10% to 90%. The 
length scale of all plots is that of the core radius of the 
potential, so that Rc = 1. 

Solid and dotted circles on the contour plots of a wave 
pattern mark the radii i?cR and i?oLR of circular orbits 
in co-rotation resonance (GR) and outer Lindblad res¬ 
onance (OLR) respectively with a neutral s = 0 mode 
with the pattern speed Dp of that mode. All pattern 
speeds are too large for there to be an inner Lindblad 
resonance (ILR). Orbits of any shape, not just circular 
ones, may be resonant. For example, the orbits which 
are in a CR with a pattern speed Dp are those which lie 
on the horizontal line D^ = SU which cuts through the 
lens-shaped region of Figure ^ from the circular orbit 
boundary on its left to the radial orbit boundary on its 
right. These orbits are spread out in space and not con¬ 
fined to one specific circle. They are more concentrated 
near that circle when most orbits are near-circular, but 
not otherwise. Similarly the orbits in an OLR are those 
for which D^ = Dp — Dij/2, and lie on a line through 
Figure n with slope —1/2. They too range from circular 
to radial. They have lower orbital frequencies than those 
in the GR because they lie further out in the disk. An 
ILR occurs only for the limited range of Dp values for 
which the line D^ = Dp -|- Dij/2 of slope 1/2 intersects 
the lens-shaped region of Figure Q] When this happens, 
there are two circular orbits, as well as a generally wide 
range of intermediate orbits, in ILR. For the most part 
we find unstable modes with growth rates s > 0. They 
have no resonances, only near-resonances at which the 
denominators of the matrix components lai are small 
when s is small. 

Below each contour plot is a plot of the radial varia¬ 
tions of the amplitude P{R) of the perturbed density (full 
curve) and of the unperturbed density (dashed curve). 
Below this is a bar chart which displays the values of the 
different Fourier components ^2 i ^^.d Wj ^ com¬ 
puted for them. The values of all these components de¬ 
pend on how the eigenvector c is normalized. We normal¬ 
ize c so as to make the positive and negative components 
of e“^®‘Dp £2 sum to ±1 respectively. We find that /C 2 ,i 
is always positive. This does not necessarily imply that 
all modes release gravitational energy and convert it to 
kinetic energy. Such a release occurs only if the sum of 
the two components IC 2 ,i + K. 2,2 is positive. Our tables 
for Kuzmin and isochrone disks list values of the ratio 
1^2,211^2,1 for each mode. A mode converts gravitational 
energy to kinetic energy if this ratio exceeds -1, and vice 
versa if the ratio is less than -1. There is no conversion 
if the ratio is exactly -1. 
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Table 1. Eigenvalues for m = 2 modes of unidirectional Miyamoto 

MODELS FOR KuZMIN’S DISK. 


riM 

-^0 


mode 



Pull Model 



/ = — 1 only 

r2p 

s 

/C2,2/^112,1 

Rcr 

l?OLR 

Qp 

3 

0 

1.000 

1 

0.825 

0.939 

1.58 

0.541 

1.296 

0.649 

3 

0 

1.000 

2 

0.418 

0.265 

1.88 

1.483 

2.246 

0.270 

5 

0 

1.000 

1 

0.913 

1.216 

1.61 

0.359 

1.176 

0.738 

5 

0 

1.000 

2 

0.530 

0.409 

1.86 

1.154 

1.878 

0.323 

7 

0 

1.000 

1 

0.963 

1.465 

2.58 

0.227 

1.115 

0.805 

7 

0 

1.000 

2 

0.643 

0.588 

2.64 

0.895 

1.609 

0.372 

3 

0.2 

0.868 

1 

1.023 

0.114 

-3.26 


1.045 


3 

0.2 

0.868 

2 

0.336 

0.203 

-1.67 

1.811 

2.632 

0.259 

5 

0.2 

0.882 

1 

1.049 

0.222 

0.53 


1.017 


5 

0.2 

0.882 

2 

0.384 

0.259 

-0.64 

1.607 

2.390 

0.294 

7 

0.2 

0.892 

1 

1.067 

0.311 

1.27 


0.997 


7 

0.2 

0.892 

2 

0.430 

0.302 

-0.17 

1.443 

2.200 

0.321 


5.1. Modes of Kuzmin Disks 

Miyamoto (1971) models are characterized by the sin¬ 
gle parameter nu- The orbital population becomes more 
nearly circular with increasing ttm, and ultimately cold 
in the limit nu oo. We follow (Hunter 1992) in work¬ 
ing with models for which all orbits circulate in the same 
direction. We use units in which G = M = Rc = 1. 

The frequencies of the fundamental modes for the self- 
consistent {Lq = 0) models listed in Tabled differ sub¬ 
stantially from those given earlier by Hunter (1992) and 
Pichon & Cannon (1997). Those results are incorrect 
because they omit the boundary integral M®. The dif¬ 
ference is considerable. Neglecting M® gives a pattern 
speed of Up = 0.357 with Rcr — 1.717 and a growth rate 
of s = 0.295 for nyi = 3. The true fundamental mode 
of the um = 3 model is the compact and rapidly grow¬ 
ing central bar shown in Figure much more compact 
than that shown in Figure 11 of Pichon & Cannon. The 
amplitude P{R) of its perturbed density in FigureEt has 
a single peak. The secondary mode, shown in the right 
panels of Figure El is slower growing and slower propa¬ 
gating. It is more extensive, has a double-peaked am¬ 
plitude, and a more spiral structure which is also largely 
confined within the CR circle. The growth rates and pat¬ 
tern speeds increase with increasing ttm- They become 
increasingly centrally concentrated as they are largely 
confined within the CR circle. 

The bar charts in Figures Et and EF display a standard 
pattern which is common to all but two of those we show. 
Angular momentum is lost by the I < 0 Fourier compo¬ 
nents, primarily ^ = — 1, and gained by the I > 0 com¬ 
ponents. All Fourier components lose 14^2,1 gravitational 
energy, much from the I = —1 component, while all com¬ 
ponents gain /C 2 ,i kinetic energy except for I = —1. The 
positive values for the ratio IC2,2/R-2,i for these modes 
show that the VV 2,2 and IC2.2 terms reinforce this trans¬ 
fer from gravitational to kinetic energy. We show in El 
how the forms of the bar charts can be understood in 
terms of the formulae we derived in a Every bar chart 
shows that a few Fourier components are significant for 


most of the transfer of angular momentum and energy. 

The second block of results in Tableware for unidirec¬ 
tional Miyamoto models from which low angular momen¬ 
tum orbits have been removed by applying the cutout 
function ll^ with Lq = 0.2. This removal reduces the 
active surface density Fact so that it tends to zero in the 
center, as shown in the dashed curve in Figure CF) but 
does not introduce a sharp barrier. It reduces the ac¬ 
tive mass Mact of the disk by a little more than 10%, 
but has a much larger effect on the modes, and espe¬ 
cially the fundamental mode. The left panels of Figured 
show the large changes in the fundamental mode of the 
nyi = 3 Miyamoto model caused by the Lq = 0.2 cutout. 
The amplitude of the fundamental mode is again single 
peaked, though its peak is moved out to i? « 0.4. This 
mode rotates so rapidly that no orbits are in CR, and it 
now extends out to the OLR circle. Its bar chart is to¬ 
tally different from that of FigureEFi but similar to that 
of Figure ini g. which is for another fundamental mode 
of a cutout disk. The I = —I components are small, 
and the flow of angular momentum and /C 2 ,i is from the 
I < 1 components, primarily I = 0, and to the I > 1 
components, primarily / = 1. 

The cutting-out of low angular momentum orbits and 
some 10% of the mass reduces the growth rate of the 
fundamental mode by such a large factor that it is no 
longer the fastest growing mode. The secondary mode, 
which is now the fastest growing, is not shown because it 
is changed much less. It has a similar but smoother spiral 
shape than that of FigureEFj and extends just beyond the 
now-larger CR circle. The two humps of its amplitude 
P{R) are displaced outwards from those of Figure EF, 
and the inner hump, which lies in a region of large mass 
reduction, is diminished. The outer hump, now at i? « 
0.9, lies in a region where the mass reduction is smaller. 
Hence the reason why the cutout affects the fundamental 
mode so much more than the secondary mode is that it 
has a much larger effect on the more central orbits which 
are the major participants in the fundamental mode than 
it does on those of the secondary mode. Because the 
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Fig. 6.— The fundamental mode (left panels) and secondary mode (right panels) of the self-consistent Lq = 0, tim = 3, Miyamoto 
model for Kuzmin’s disk; the first two entries of Table [Tl Here and in all subsequent plots of modes, the top panels show positive contours 
of the perturbed density Si, in steps of O.lSi. The solid and dotted circles mark the co-rotation and outer Lindblad resonance circles. 
The middle panels show the wave amplitude (solid line and left scale) and unperturbed density (dashed line and right scale). The bottom 
panels show Fourier components of kinetic energy (grey bars), gravitational energy (white bars), and angular momentum (thin bars). 


ratio K, 2 ^ 2 /K, 2 ,i < —1 for the both modes of the cutout 
riM = 3 model, both induce a transfer from kinetic to 
gravitational energy. 

The results for higher nu are similar. The cutout 
increases the pattern speeds of fundamental modes so 
much that none have CRs. It decreases the pattern 
speeds of secondary modes, though not enough for there 
to be ILRs, which are possible only when flp < 0.130. 
All growth rates are decreased, though the fundamen¬ 
tal mode is still, though barely, the faster growing for 
Lq = 0.2 and nu = 7. As nu increases and orbits be¬ 
come more circular, the relative amount of angular mo¬ 
mentum absorbed by the I = 1 component of secondary 
modes increases. 

Toomre’s axisymmetric stability parameter Q in¬ 
creases monotonically outwards from its central value of 
27r/[3.36(l -I- nM/2)^/^] for Miyamoto models. Kalnajs 
(1976b) introduced an alternative set of models, which 
also depend on an integer parameter tok which increases 
as the models cool, and for which Toomre’s Q remains 


almost constant. Sellwood & Athanassoula (1986) modi¬ 
fied these models with the addition of two extra parame¬ 
ters: an angular momentum Jc and /3. They reversed the 
sense of rotation of a fractional mass Mretro of stars in a 
unidirectional Kalnajs model with angular momenta in 
the range (0, Jc), following equation (5) of Zang & Hohl 
(1978). This gives a smoothly tapered DF, and elimi¬ 
nates the discontinuity of the unidirectional model. Our 
results for some of these models are given in table |21 and 
displayed in Figures 0 and |H1 The fundamental mode of 
the unidirectional ttik = 6 model, listed in Table El has 
not been plotted because it is a compact and rapidly ro¬ 
tating bar, a little larger, but otherwise just like that of 
the Miyamoto model shown in the left panels of FigureEl 
The right panels of Figure 0 show how the compact bar 
is modified by a tapering with Jc = 0.25 which reverses 
less than 5% of the orbits. Like the cutout, it gives the 
fundamental mode a smoother, more spiral, and more 
extensive pattern, and a slower growth rate. A major 
difference is that the taper diminishes the pattern speed. 
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Table 2. Eigenvalues for m = 2 modes of Athanassoula & Sellwood 

MODIFIED KALNAJS MODELS FOR KuZMIN’S DISK 


ruK 

P 

Jc 

-^retro 

mode 


Full Model 


AS 

a 

P04*> 

/ = — 1 only 


S K-2,2 

IK.2,1 

Op 

S 

Op 

s 

Op 

4 

0 

0.40 

0.084 

1 

0.335 

0.174 

1.83 

0.168 

0.020 



0.193 

6 

0 

0 

0 

1 

0.746 

0.711 

1.42 





0.569 

6 

0 

0.25 

0.046 

1 

0.445 

0.308 

0.52 

0.233 

0.066 

0.22 

0.114 

0.264 

6 

0 

0.25 

0.046 

2 

0.294 

0.109 

1.85 

0.165 

0.058 

0.175 

0.055 

0.207 

6 

3 

0.60 

0.154 

2 

0.158 

0.027 

-4.00 

0.145 

0.014 

0.14 

0.02 

0.144 

8 

4 

0.90 

0.160 

1 

0.199 

0.064 

-1.47 

0.173 

0.035 



0.174 


^Athanassoula &; Sellwood (1986) 
’^Polyachenko (2004) 







Fig. 7.— The left panels are for the fundamental mode of a cutout riM = 3, Lq = 0.2, unidirectional Miyamoto model. The right panels 
are for the fundamental mode of a tapered mx = 6, Jc = 0.25, Athanassoula & Sellwood model for which some 4.6% of the orbits have 
been made retrograde. 


whereas the cutout increases it. Both the tapering and 
the cutout cause a large change in the population of low 
angular momentum orbits which dominate in the central 
part of the disk where the fundamental mode is concen¬ 
trated. Whereas the cutout removes many of them, the 


tapering merely reverses the rotation of many, and makes 
the DF isotropic for \L\ <C Jc- Figure [7|f shows a bar 
chart which is of the standard pattern, though the ta¬ 
pering has increased the magnitudes of the I > 0 energy 
components. 
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Fig. 8.— The left panels show the fundamental mode of the rriK = 4 Athanassoula &; Sellwood model, tapered with Jc = 0.4 and with 
8.4% of its orbits retrograde. The right panels show our computed mode for their /3 = 3, mx = 6, Jc = 0.6 model with 15.4% retrograde 
orbits. 


We computed the fundamental mode for the tapered 
rriK = A, Jc = 0.4 model because it is one of the few 
for which Sellwood & Athanassoula plot a mode shape. 
Our results in the left panels of Figure |H1 are similar to, 
though less spiral than, the fundamental mode of the 
cooler and more sharply tapered disk in the right panels 
of FigureO though there is a remarkable dearth of angu¬ 
lar momentum transferred to its Z = 1 component. The 
major difference between their Figure 5b and our Figure 
IHh is that theirs, which is less spiral, has some central 
structure which is not present in Figure |Hh- 
Table El compares our results with those of Athanas¬ 
soula & Sellwood (1986) for four of the 33 models for 
which they did A^-body simulations. Some discrepancies 
between the matrix theory and the simulations with their 
finite number of particles and gravity softening are to be 
expected. Nevertheless, the large discrepancies between 
our results and theirs, and that of Polyachenko (2004), 
for the (3 = 0 models are worrisome. Our best explana¬ 
tion for it is that we have found these disks to be very 
sensitive to near-radial orbits, and, as we explain in an 
and ^4.51 we find eigenvalues to be sensitive to the accu¬ 


racy with which the central regions are handled. 

In contrast to these discrepant cases, our results for the 
two non-zero (3 models listed in TableElagree reasonably 
well with those of Athanassoula & Sellwood. A positive 
value of (3 removes orbits with low energies and angular 
momentum high relative to Lc{E), and adds to those 
with higher energies and lower angular momentum. It is 
easy to show from their Appendix A that 


fo{E,0) = 


{mK-2(3El) ^ RcE 


27r^ 


■exp[(/3(l-F,^)], Es = 


GM ■ 
(78) 

The value of (3 can not exceed m]^l2 because /o is then 
negative at the center of the disk where the scaled en¬ 
ergy Eg = —I. Many of Athanassoula & Sellwood’s mod¬ 
els have marginal values (3 = —rriKl^. Their /o values 
peak at intermediate values of Eg, and they have dou¬ 
ble peaked radial velocity profiles. The {rriK, (3, Jc) = 
(6,3,0.6) model is one of two to which Polyachenko 
(2004) applied his simplified theory to get a value of the 
pattern speed Qp which agrees with what we get from 
I = —1 terms only. Our plot of the most unstable mode 
for this case in Figure Et shows the same concentrated 
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central structure as in Polyachenko’s Figure 6, but his 
plot has a stronger second hump around R — 1 than our 
weaker one, and lacks the extended outer spiral that we 
find. The structure of this mode suggests that it is a 
secondary one, though it is the only one we have found. 
Its low pattern speed means that it comes closer to hav¬ 
ing an ILR than any other mode in our survey. Its bar 
chart in Figure IHF is noteworthy for its large and nega¬ 
tive I = —1 component of /C 2 p. This mode is another one 
for which K. 2 . 2 /K. 2 ,i < —1, and hence converts kinetic to 
gravitational energy. 


5.2. Modes of Isochrone Disks 

Kalnajs (1976b) gives unidirectional models for the 
isochrone disk. They contain an integer parameter tok 
which increases as the models cool, and have fairly uni¬ 
form values of the Toomre parameter Q. Kalnajs cal¬ 
culated modes for modified versions of these models in 
which the sense of rotation of some stars is reversed to 
make them retrograde. Kalnajs’s formula for retrograde 
stars is given in equation (13) of Earn & Sellwood (1995). 
Combining that with the x —> 0 limit of equation (26) of 
Kalnajs (1976b) and using GR formula (7.126.1) gives 


/o 


_TOK_ 

67r2 [1 + J^ + |J^|]2 -k-2- 


< 0 , 


(79) 


for the DF of retrograde stars in units in which G = M = 
Rc = 1. There are now retrograde stars of all angular 
momenta, not the limited ranges of the models of EH 
Integration over phase space gives 


A/petro — 


_Tok_ 

3(tok — 2)(2mK — 3) ’ 


(80) 


for the fractional mass in retrograde stars. 

TableElgives our results and those of others in units in 
which G = M = Rc = 1 for two modes of these models. 
The mutual agreement is now gratifyingly close. Both 
pattern speeds and growth rates increase with increasing 
mK, as they do with Miyamoto models with increasing 
riM- Pichon & Cannon (1997) found secondary modes 
for some values of tuk other than those we have listed 
in our table. Again there are no ILRs. ILRs occur only 
when Up < 0.0593 in the units used here, for which the 
ranges of both and are only a half of what they 
are in Figure ^ 

Figure El displays two modes of the isochrone disk with 
Tok = 12. The wave patterns match those of Earn & 
Sellwood’s Figure 1. Spirality increases with tok- The 
secondary mode is more extended than the fundamental, 
and the amplitude of its spiral arm is two-peaked, ver¬ 
sus the single peak of the fundamental mode. Both bar 
charts show that the I = 1 component absorbs much more 
angular momentum and kinetic energy than the I = 0. 
This is different from Figure El but is in line with the 
trend we noted with Miyamoto models with increasing 
riM, though note that tuk = 12 corresponds to a much 
larger um value than the riu = 3 of Figure El The sec¬ 
ondary mode of the warmest tuk = 6 model is the only 
one for which kinetic is converted to gravitational energy. 

We studied some other models for the isochrone po¬ 
tential, and found, as did Pichon & Cannon (1997), 
that their modes are qualitatively similar to those for 
Kuzmin’s disk. Only their scales differ because the 
isochrone disk is the more spread out. 


5.3. Modes of the Exponential Disk 

Our exponential disk models are the most varied. After 
normalizing with units in which G = vq = Rc = 1, we 
are still left with three parameters; N which measures 
the tendency of the orbits to circularity (cf Figure El) 
i?D which measures the length scale of the exponential 
disk relative to the core radius of the potential, and the 
density scale E^. A cutout introduces a further parame¬ 
ter Lq. The central density of the disk is E;,e“''' and its 
total mass Md = 27rEsi?D(7?c + i?D)e“'^. As FigureEh 
shows, disks become increasingly centrifugally supported 
with increasing N and hence this parameter is similar to 
the parameters um and mK which we varied for previous 
models. Here we restrict our exponential disk models to 
the single case of A^ = 6 to allow us to study the conse¬ 
quences of varying the other parameters. As Figure Et 
shows, N = 6 models have Toomre’s Q parameter close 
to unity over a substantial central part of the disk, whose 
size increases with Ro when Rd > 1. 

TableEJlists the fundamental and secondary modes for 
six different N = 6 models. It lacks values of the ratio 
1 ^ 2 , 2 / 1 ^ 2,1 because we were unable to compute accurate 
values of JC 2,2 = —kV 2,2 wit h ou r chosen basis functions, 
for the reasons discussed in EH ILRs occur in the cored 
logarithmic potential only for Up < 0.106, which is much 
smaller than any of the pattern speeds. Figure 11171 shows 
the two modes for the first Rc = 1, S^Rd = 0.42, model. 
The fundamental mode is a rapidly rotating bar confined 
within the CR circle, and similar in all respects to that of 
the Kuzmin disk in Figure Eh- The secondary mode has 
a slightly lower pattern speed and growth rate, and has a 
spiral form. Its amplitude has the usual two peaks within 
the CR circle, and its pattern extends a little beyond the 
CR circle but lies within the OLR circle. Its bar chart is 
remarkable for its small I = 0 components. 

The first four models of Table 0 are chosen to study 
the effect of varying Rc- The parameter HsRu has to be 
changed too because it is necessary to remain within the 
physically allowed region of Figure El Its values are near 
marginal in that they are approximately 90% of their al¬ 
lowed maximum. Despite the decrease of 'E^Rc, the total 
disk mass grows as its scale length Rd increases, though, 
as Figure El shows, the halo/bulge grows in importance 
as i?D decreases, and the disk becomes progressively less 
maximal. The structure of the fundamental bar mode 
does not change along this sequence; it remains a com¬ 
pact and rapidly rotating central bar. Figure lllb for 
Rd = 1-6 shows that the amplitude of the secondary 
mode has developed a third hump. This development 
occurs around Rd « 1.25, and seems to be related to the 
anomalously low growth rate at Rd = 1-2 (See TableEJ. 
The biggest difference between the bar chart Figure [TTV: 
for the larger Rd = 1.6 disk and that of Figure If OH for 
Rd = 1 is the greatly increased significance of the I = 0 
component. 

Figures Ea and ESI illustrate how decreasing Eg and 
hence the mass of the disk, whilst keeping its length scale 
Rd fixed, stabilizes the disk. The transition of the fun¬ 
damental mode from stability to instability appears to 
take place through a pitchfork bifurcation. That of the 
secondary mode appears to take place through a tan¬ 
gent bifurcation. The order in which the two modes are 
stabilized is different for the two different values of Rc- 
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Table 3. Eigenvalues of Kalnajs models for the isochrone disk. 


m-K 

-^retro 

mode 


Full Model 


K78^ 


/ = — 1 only 

Qp 

S K.2,2 

//C2,l 

Qp 

s 

f^p 

s 

f^p 

6 

0.056 

1 

0.169 

0.080 

-0.11 

0.170 

0.075 

0.170 

0.075 

0.094 

6 

0.056 

2 

0.121 

0.035 

-1.58 





0.079 

9 

0.029 

1 

0.235 

0.149 

-0.05 

0.235 

0.145 

0.235 

0.145 

0.127 

9 

0.029 

2 

0.183 

0.089 

-0.92 





0.109 

12 

0.019 

1 

0.292 

0.217 

0.07 

0.295 

0.210 

0.295 

0.210 

0.150 

12 

0.019 

2 

0.234 

0.148 

-0.60 



0.230 

0.145 

0.133 


^Kalnajs (1978). 

’^Pichon & Cannon (1997). 



When Esi?D is decreased at a fixed Rd, then the part of 
the rotational velocity due to the halo/bulge components 
increase from those shown in Figure 13 which are drawn 
for values of T,sRb which are 90% of the allowed max¬ 
imum. The stabilization shown in figures 1121 and na is 
therefore similar to that which is achieved by sufficiently 
massive halos (Kalnajs 1972, Ostriker & Peebles 1973, 


Hohl 1975). 

The stabilization of both modes in the neighborhood 
of TisRd ~ 0.3 for both Rd = 1 and Rd = 1.6 suggests 
that the stability boundary approaches the boundary of 
the physically feasible region plotted in Figure El as Rd 
increases, i.e. as A decreases. It raises the possibility 
that the two boundaries intersect before the A ^ 0 limit 
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Table 4. Eigenvalues for m = 2 modes of exponential disks with N = 6. 


mode 

Ro 


-^0 

-^act 


Full Model 


[ = — 1 only 


s 

^CR 

^OLR 


1 

1 

0.42 

0 

1.000 

0.914 

1.151 

0.444 

1.688 

0.698 

2 

1 

0.42 

0 

1.000 

0.693 

0.324 

1.040 

2.334 

0.364 

1 

1.2 

0.38 

0 

1.000 

0.840 

0.935 

0.646 

1.871 

0.591 

2 

1.2 

0.38 

0 

1.000 

0.607 

0.059 

1.309 

2.701 

0.330 

1 

1.4 

0.36 

0 

1.000 

0.805 

0.794 

0.737 

1.967 

0.540 

2 

1.4 

0.36 

0 

1.000 

0.467 

0.179 

1.893 

3.572 

0.294 

1 

1.6 

0.34 

0 

1.000 

0.768 

0.642 

0.834 

2.077 

0.488 

2 

1.6 

0.34 

0 

1.000 

0.443 

0.119 

2.024 

3.775 

0.277 

1 

1 

0.42 

0.1 

0.967 

1.170 

0.259 


1.212 


2 

1 

0.42 

0.1 

0.967 

0.537 

0.288 

1.571 

3.082 

0.331 

1 

1 

0.42 

0.3 

0.913 

1.044 

0.228 


1.423 


2 

1 

0.42 

0.3 

0.913 

0.439 

0.274 

2.047 

3.811 

0.302 







Fig. 10.— Modes of the exponential disk for N = 6, Ry) = 1, and SsKd = 0.42. The left panels are for the fundamental mode, and the 
right panels are for the secondary mode. 
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Fig. 11.— The secondary mode of an Af = 6 exponential disk of 
larger extent with Rj^ = 1.6 and SsRu = 0.34. 



Fig. 12.— The variation of the pattern speed flp (open squares) 
and growth rate s (filled squares) of the N = 6, Ry) = 1 exponential 
disk as the mass of the disk is varied. Data for the fundamental 
modes (1) and secondary modes (2) are connected by solid and 
dashed lines, respectively. The secondary modes are double peaked. 



Fig. 13.— Same as Figure IT^ but for N = 6, Ry, = 1.6. The 
secondary modes now have triple peaks. 


of an exponential disk in a singular logarithmic poten¬ 
tial, for which SsRd = 0.304 is reached. Our current 
computer algorithms are not capable of approaching that 
limit, but they do show the Ru = 2 disk to be quite 
stable. This suggests that the classical exponential disk 
(A = 0) with a completely flat rotation curve (i?c = 0) 
may be stable against bisymmetric excitations. This disk 
is much less than maximal, and requires a substantial 
central bulge/halo to maintain its equilibrium. 

The final Figure IT^ and the last two models of Table 
^show the effect of cutting out low angular momentum 
orbits via the cutout function There is a major 

change between Lq = 0 and Lq = O.I. Comparison with 
Figure irni shows that the fundamental mode is changed 
much more than the secondary one, as we found with the 
Kuzmin disk. Its growth rate is diminished substantially, 
and its pattern speed is increased so much that there is 
no longer a CR circle. It is more spiral and extensive 
as the peak of its amplitude has moved outwards from 
the region in which the density has been diminished sub¬ 
stantially. Its bar chart has undergone a large change and 
now resembles that of Figure |7^; the / = — I components 
are insignificant and the flow of both angular momentum 
and K 2.1 is primarily from Z = 0 to ( = 1. The growth 
rate and pattern speed of the more extensive secondary 
mode is changed much less by the cutout. The decrease 
in the inner density of the disk leads to the near total 
disappearance of the inner hump of its amplitude. The 

1 — 0 components are again significant in the bar chart of 
Figure ITlif. unlike the i?D = 1 case of Figure UnK but like 
that for the Rd = 1-6 case of FigureITTfc. The modes for 
Lq = 0.3 are similar to those for Lq = O.I. The singular 
behavior of e~ ' ° makes the study of small values of Lq 
and the approach to the limit Lq 0 computationally 
difficult. 

We can not compare our results directly with those 
of Sellwood (1989) which are for uncored exponential 
disks in the cored logarithmic potential lEOJl. Not only 
are the disks different, but his are for the larger range 

2 < Ryq/Rq < 8.33 of disk radii than ours. His m = 2 

modes have pattern speeds that are also too fast for ILRs, 
but the spiral mode for Rd/Rc = 5 shown in the left of 
his Figure 3 extends out to the OLR circle, and so is more 
extensive than any of ours for exponential disks without 
cutouts. As noted in EH Sellwood’s value of 0.6 for the 
critical parameter uo(Rd/G'Md)^/^, which corresponds 
to [e^/27rSsRD(l-b in our scalings, is smaller than 

any of ours. That critical parameter is 0.86 when sta¬ 
bility is achieved for EsRd = 0.29 in Figure IT^ for ex¬ 
ample. That value is still significantly less than the 1.1 
which Efstathiou, Lake & Negroponte (1982) found to be 
necessary for their A^-body experiments, but that larger 
value may well be necessary for nonlinear stability. 

6. DISCUSSION 

We analyze and explain the properties of the bar charts 
in EH and discuss the implications of our findings for 
Polyachenko’s (2004) simplified theory in ^16.21 

6.1. Transfer of Angular Momentum and Energy 

The rate of change of the perturbed angular momen¬ 
tum (IA 6 II contains a factor s/\lLl}i + mVl,j, — tupp, which 
tends to tt [m(D 0 — Lip) -b ILIr] as s ^ 0. As Lynden-Bell 
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Fig. 14.— Modes of the TV = 6, /2 d = 1? and Ssi?D = 0-42 exponential disk with an Lq = 0.1 cutout. The active surface density is 
shown by the dashed lines in the central panels. The left panels show the fundamental mode and the right panels the secondary mode. 


& Kalnajs (1972) note, this implies that a neutrally sta¬ 
ble wave emits and absorbs angular momentum only at 
resonances. For the unstable modes we discuss, the fac¬ 
tor s/\inR -I- — wpp peaks at resonances, sharply 

so when s is small. This peaking implies that emission 
and absorption occurs mainly near resonances. Figure ^ 
shows that, for given pattern speed flp > 0, the orbits as¬ 
sociated with I > 0 resonances lie at successively smaller 
values of the orbital frequencies flu and and hence 
successively further out in the disk. That is the reason 
for interpreting the flows of angular momentum, and also 
kinetic energy, which are proportional to the quantities 
shown in the bar charts, as outward flows. 

It is evident from equations 123), l|S3i and that 
the partial derivatives dfo/dJu and dfo/dJ^ play a ma¬ 
jor role in determining the signs of the Fourier compo¬ 
nents shown in the bar charts. All the other terms in 
7^2 are either magnitudes or constants. Lynden-Bell & 
Kalnajs argued that both derivatives are negative for 
a physically reasonable DF. If this is so, then equation 
shows that L 2 > 0 toi I > 0. The compensating 
negative values of necessary because there is no net 
change of angular momentum, can and must occur for 


I < 0. The fact that is negative, as we find, can 
be accounted for by \dfo/dJii\ being much larger than 
\dfo/dJ^\. The prediction matches the standard pat¬ 
tern we find in The assumptions that dfo/dJa < 0 
and \dfo/dJR\ \dfo/dJ^\ are generally valid. Fig¬ 
ure 6 of Kalnajs (1976b) gives examples; the contours of 
the two unidirectional DFs shown there decrease much 
more rapidly with increasing than with J^. However 
the modified /3 = mKl‘2 models of H5.ll are exceptional 
because they have regions in which dfQ{E,L)/dE, and 
hence dfo{JR, J^)/dJR, are positive. 

Though positive values of dfo/dJn are unusual, pos¬ 
itive values of dfo/dJ^ are not. They are unavoidable 
with disks which are totally, or mostly, unidirectional. If 
fo {Jri 0) = Oj it is for the unidirectional models of 
Zang (1976) and Evans & Read (1998a) and our cutout 
disks, then dfo/dJ^ must be positive for some positive 
values of because otherwise could never 

become positive for > 0. Similarly dfo/dJ^ must be 
positive for some > —Jc for the tapered models of 
Table H Furthermore the explicit formula shows 
that dfo/dJ^ > 0 ior < 0 for the isochrone models of 
H5.2I Although no positive values of dfo/dJ^ are visible 
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in Figure 6 of Kalnajs (1976b), these are disks for which 
the boundary integral 


2 2 
— m TT 


dJ_R 




J<p—0 


(81) 


must then be included in the expression for It is 
negative for every I, as are the contributions to the in¬ 
tegral m from positive values of dfo/dJ^j). Unidirec¬ 
tional disks for which Jq {Jr, 0) > 0 can be regarded as 
Jc ^ 0 limits of tapered disks; the boundary integral ED 
then accounts for the effect of the large positive values 
of dfo/dJ^ which occur in the taper. 

Despite the occurrence of regions of positive values of 
df[)/dJ(f,, they are generally either too small, or confined 
to too limited regions, to modify the standard bar chart 
pattern. However, cutout disks have central regions in 
which dfo/dJ^ is large and positive. That is the reason 
for the negative values of L® f^e barcharts of Figures 
and [T^ for fundamental modes of cutout disks. The 
bar charts for the secondary modes of these cutout disks 
have the standard form, as is seen in Figure EK This is 
because dfo/dJ^ < 0 for the more distant orbits which 
are the more important for the secondary mode. 

Integral 18811 for i differs from integral H.'ITII for L 2 
by a factor (D^ -|- IVIr/2). This is positive for I > —1 
for all direct orbits, but negative for / < — 2 (See Figure 
nj, and is the reason why the signs of 1 match those 
of L 2 for ^ ^ but are their opposites for I < —2. 
This pattern holds even for the exceptional cases, with 
the result that AT^ 1 is positive for all I, except for I = —1 
in standard bar charts and Z = 0 in our exceptional ones. 

Integral 18511 for IF 2 1 differs from that for L 2 by a fac¬ 
tor (Dp — — IVIr/2). This factor is positive for all 

orbits for I = 0 modes when the pattern speed Dp > I, 
which is why W 2 1 is negative like L 2 for the exceptional 
fundamental modes. IU 2 1 s-^rl A® otherwise have oppo¬ 
site signs for other modes with slower pattern speeds and 
A 2 > 0, and which are dominated by orbits within CR. 
In fact IU 2 I is negative for all I in all bar charts. This is 
explained by noting that (Dp — D^ — IQ.r/2) is negative 
for modes which lie principally within OLR for Z > I, 
and positive for Z < —I, so that the signs of ^ are 
respectively the opposite and the same as those of A^- 

A striking feature of the bar charts is the rapidity with 
which the quantities displayed there decrease with in¬ 
creasing |Z|. There are two reasons for this. One is the 
increase in the denominators of their integrands in the 
regions where most orbits lie. The Z = — 1 terms are gen¬ 
erally important, even though there is no ILR, because 
the denominator term 2(D,^ — Dp) — D/j is never large 
for any orbit. Another reason is the decay of the Fourier 
coefficients Vi with increasing |Z|. 


6.2. Abnormal Orbits and Polyachenko’s Theory 

Tables H through 0] show consistently that eigenval¬ 
ues estimated using Polyachenko’s simplified theory are 
lower than those calculated from the full matrix the¬ 
ory. Pattern speeds Dp well exceed the narrow range 
of Di values which the cored potentials considered here 
allow, and so his basic assumption of small |Dp — D^j is 
then not valid. Our bar charts show clearly that a few 


Fourier components other than Z = — 1 are always signifi¬ 
cant. The orbits which participate in the unstable modes 
that we find are predominantly abnormal in the sense de¬ 
fined by Lynden-Bell (1979) because all sufficiently cen¬ 
tral, as well as more radial, orbits are abnormal. In fact 
even circular orbits are normal only for R/Rc > 2.44 
for Kuzmin’s disk, R/Rc > 3.73 for the isochrone, and 
R/Rc > 25.3 for the cored logarithmic potential (See 
Figure 1^. Our modes lie primarily within these limits. 

A quite different situation arises with the scale-free 
potentials Vo{R) = sgn{a)R°‘ for a ^ 0, and Vo{R) = 
Ini? for a = 0, studied by Evans & Read (1998a,b). 
Even central orbits can then be normal because the 
normal/abnormal classification then depends solely on 
the ratio y = L/Lc{E), and is independent of energy. 
Lynden-Bell’s criterion lEll is satisfied, and an orbit ab¬ 
normal, if 


fiy) d 

'2-a 


7/y 

2 


fiy) - 


a 


ny) + ^\ [fiy)-yfiy)]>0, (82) 


where the function f{y) = ttJr/L c{E) and depends also 
on a. Its derivative f'{y) = —ttD^/D/j = ^\g{a,y)\, 
where g{a,y) is the function defined and analyzed in 
Touma & Tremaine (1997). The criterion H8‘ill gives 
y < 0.723 for the scalefree logarithmic potential, so that 
all the more circular orbits with large L/Lc{E) are nor¬ 
mal. For a = —0.25 and a falling rotation curve, there 
is a wide range of normal orbits, and only orbits with 
y < 0.496 are abnormal. For a = 0.25 and a rising ro¬ 
tation curve, only nearly circular orbits are normal, and 
all those with y < 0.973 are abnormal. All orbits are 
abnormal for a > 0.275. Interestingly Evans & Read, 
who looked specifically at the cases of a = ±0.25, found 
that unstable modes grow more vigorously for the rising 
rotation curve case of a = 0.25 with many abnormal or¬ 
bits. ILRs occur with their modes because the range of 
Di is unbounded for scale-free potentials. 


7. SUMMARY 

This paper develops the theory of modes in thin stellar 
disks. It then implements that theory for a selection of 
disks. The theory basically is that pioneered by Kalnajs 
(1971,1977). Our development of its Eulerian form to 
second order is new. The expressions for angular mo¬ 
mentum and total energy were given earlier by Kalnajs 
(1971) and Lynden-Bell & Kalnajs (1972), but the ex¬ 
pressions for the two separate components are new. The 
boundary integral terms which must be included in the 
method are also new. These integrals are avoided when 
a Lagrangian form of the theory is used. The Lagrangian 
form leads to more complicated computations, which is 
why prior computational work, except for that of Van- 
dervoort (1999), has all been done for using the Eulerian 
form. We show in Appendix H how the passage from 
the Lagrangian to the Eulerian form allows us to inter¬ 
pret the boundary integrals as boundary fluxes. This 
analysis does not apply to the quite different version of 
Lagrangian theory which Vandervoort uses. 

Our applications of the theory are to m = 2 modes in 
potentials with smooth cores. In all but one case (see 
Tabled, our search has found two modes, a centrally 
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concentrated fundamental mode with a single peak in 
amplitude, and a more extensive and more spiral sec¬ 
ondary mode with at least two peaks. There may be oth¬ 
ers. We find that the shape, pattern speed, and growth 
rate of the fundamental mode are especially sensitive to 
the relatively small proportion of orbits which provide 
the density in central regions because they have low an¬ 
gular momentum. The growth rate of the fundamen¬ 
tal mode is reduced substantially, and its shape becomes 
more extensive and spiral, if those orbits are either re¬ 
moved or reversed. Removing them increases the pat¬ 
tern speed, while reversing them decreases it. Removing 
them also makes the more extensive and more spiral sec¬ 
ondary mode the faster growing in all but one of the 
cases in Tables ^ and 0| This sensitivity to orbits, many 
of which are near radial, might suggest a connection with 
the phenomenon of radial orbit instability. However that 
phenomenon, reviewed recently in Merritt (1999), occurs 
in hot anisotropic spherical stellar systems in which ra¬ 
dial orbits are sufficiently predominant. The sensitivity 
we find here arises even with cool stellar systems with 
few near radial orbits. 

With two exceptions, modes are largely confined within 
the CR circle, but are too fast for there to be any or¬ 


bits in ILR. The lack on an ILR means that modes can 
propagate into, and be reflected from, the center of the 
disk, even in cases in which our densities drop to zero 
there. The modes are unstable, sometimes rapidly so, as 
swing-amplifier theory (Toomre 1981) predicts. We have 
achieved stability only by decreasing the active mass of 
the disk. In the exceptional cases, the modes lie within 
the OLR circle, and the pattern speeds are then too fast 
for any orbits to be in CR. All the unstable modes trans¬ 
fer angular momentum and kinetic energy outwards, and 
most release gravitational energy and convert it to kinetic 
energy. This flow of angular momentum and energy is de¬ 
rived from the second order extension of a linear theory, 
and so can describe only the early stages of the growth of 
an instability, and not its later fully nonlinear develop¬ 
ment. Only a few Fourier components account for almost 
all of the angular momentum and energy. Polyachenko’s 
(2004) theory, which is equivalent to ignoring all but one 
term of our Fourier development, seems to be an over¬ 
simplification. 


This work has been supported in part by the National 
Science Foundation through grant DMS-0104751. 
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APPENDIX 

ANGULAR MOMENTUM AND ENERGY 


Both the perturbed angular momentum £2 and the first component JC 2 ,i of the kinetic energy contain integrals of 
the form 

l2{t)= JJ SiJR,J^)f2d3dQ, (Al) 

for different functions S. Such integrals can be evaluated as follows 

= JJ S{Jn, J^)^d3de = JJ S{Jr, J^) + [/2,7fo]) dJdO, (A2) 

because the terms added are angle derivatives of periodic functions which integrate to zero over the angles. Equation 
m gives the term in parentheses as the sum of two terms. The first contributes 

/f = /f =JJik 


Here the subscript i/ represents the pair of subscripts {R,(p), and we suppose the summation convention to apply to 
it. The last step again uses the fact that any integral of an integrand which is a derivative with respect to an angle 
vanishes. The second term contributes 


JJ -S[h,V,]d3dQ 


dfi aui 

dJ, 80, 


dfi dVi 
80, 8J, 


dJde = 


8Vi 8S 8Vi 8{Sh) 
80, 8J, ^ 80, 8J, 


djSh) 8Vi 
80, 8J, 


dJde. 


The combination of the second and third components vanishes because it can be combined to an integral of a divergence: 



8 

W, 



dJde. 


(A5) 


The second set of component with angle derivatives integrate to zero, but so too do the integrals of derivatives with 
respect to the actions. That is because the differentiated quantities vanish at the limits in action space, as —*■ ±00 
and ^ 00 where the perturbation tends to zero, and at J/j = 0 where Vi is independent of because 4'™ (0, J^) = 0 
for Z ^ 0 (cf H2.2II . The remaining first component of (IA4II can be evaluated for S = as in El, to give 


dC2{t) 

dt 


—2ms7r^e^®* ^ dJ (l 


1— — C 


dfo 8fo \ 
-1- m - I 




8Jr dJ^J \K}r + - uj\ 


(A6) 


This result agrees with that of Lynden-Bell & Kalnajs (1972) when account is taken of the fact that their Fourier 
coefficients are larger than ours by a factor of 47r^, and their uj has the opposite sign. For S' = Tfo, we get 


^ / / Ho/2dJd0 = 


1 — — OC 


dJ I 


dfo 

8Jr 


A simple consequence is that 


d£’2 d(/C2+W2) 


dt 


dt 


= D, 



(ZUij + mD0)|V)p 


\lflR + — wp 

dC2 


' dt ’ 



(A7) 


(A8) 


where £ is the total energy. The undifferentiated values which are quoted in ES follow because of the simple time 
dependence on 

A deeper analysis of the second order equation ®, though not its full solution, is needed to evaluate the integral 


W2,2 = = JJ Vo/2dJd0. (-^9) 

We rewrite equation m as 

^ + [f2,no] + [fo,V2] = -[/i,Ui]. (AlO) 

The left hand side of irrm . which is homogeneous in subscript 2 quantities, has the same form as the first order 
problem for which we derived the homogeneous linear equations Equation Einii leads in a similar way to 

inhomogeneous linear equations. Its right hand side contains both axisymmetric terms and non-axisymmetric ones 
with wavenumber 2m. We need consider only the axisymmetric terms and the part of the solution for /2 which they 
cause, because only they will contribute to the integral 1A9II for yV 2 , 2 - They have a Fourier expansion 


00 

g 2 si ^ where 

I —— 00 


1 




(All) 
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We represent the potential and density of the axisymmetric part of /2 by series 


V 2 = ^ S 2 = ^ a,a°(E), (A12) 

j=o j=o 

like those of equations 03 and ca but now in a complete set of axisymmetric basis functions. We use Fourier 
expansions 


/2 = e'*‘ E U2 = e2*‘ Ui{JR,J^)e 


UBr 


1 — — C 


1 — — C 


like those of equations cu and with 




(A13) 


(AM) 


j=o 


and with Fourier coefficients defined as in equation m for m = 0. Matching Fourier coefficients in equation 
irm gives 


(2s + imR)gi-il^Ui = Ni. 

oJr 

Then, following the same procedure as used in EH we obtain the matrix equation 

[M(0,2iis) - D(0)]a = h, 

where the components of the column vector h are given by 

I (/O D — 

I —— 00 ‘ 


= 4^' E / 


{IflR — 2is) 


(A15) 

(A16) 

(A17) 


The matrix M(0, 2is) is real because each ±Z pair in the sum 121 II comb ines two complex conjugate quantities, be cause 
is even in I [cf eq. II 8 I 1 I. The right hand side h of equation IAI 6 II is real because the ±1 pairs in the sum 1A17II 

also combine two complex conjugate quantiti es, du e also to the fact that Ni = N-i because the Ni are the Fourier 
coefficients of a real function. Hence equation ITT^ is a real matrix equation, and its solution for the unknown vector 
a is real. Knowing a, we can evaluate 


00 

/ OO n 

VbS 2 dx = 27re2"* y] Oj / Vo{R)a^{R)RdR. 

j=o i 


(A18) 


The reason why it is so much easier to compute 14^2,1, ^^ 2 , 1 , and C 2 is that they need only the single Fourier coefficient 
qq. Equation gives go simply as Nq/2s and the matrix equation is not needed. 

The computation of a can be checked by verifying that the total mass due to the axisymmetric density F 2 vanishes. 
This is 


00 „ 

M2it) = 2-Ke^^^ Qj / a°{R)RdR = 0. 
f=o n 


(A19) 


Formally, the constancy of Af 2 follows fro m the analysis of Appendix^ it is the S' = 1 case of integral lAlll . With 
the Glutton-Brock functions lHiHll . the sum 1A19II can be evaluated as using GR formula (7.225.3). 

Substituting the Fourier series m and m for /i and Vi and carrying out the angle integrations gives 


= i E 


k— — c 


[fkVk-i - fkVk+i) - I Ivk-i^ + Vk+i 


dfk 

OJr 


(A20) 


Derivatives of fk, and hence second order derivatives of /o, are avoided by integrating equation iXlTll by parts with 
respect to the actions. This gives the following integral over the whole action space: 
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For the unidirectional disk with DF (1^ . hj is given by the integral (IA21II over > 0, > 0 with /o = /(f, plus 

the following boundary integral: 
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where 
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(A23) 


kQfi + mVL^ — uj kVLR + — ui 

Note that the solution for yV 2,2 = —V. 2,2 intermingle different Fourier components, unlike £ 2 , W 2 ,i, and /C 2 ,i for which 

the Fourier components can be separated as in equations and 123- _ 

The partial derivatives of T™ (J/j, J^) needed for equations (IA21II and IIA22II can most easily be calculated simulta¬ 
neously with the '^Y^AJr, J^). For this we differentiate equation 11811 partially with respect to an action to obtain 
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The variables i?, (0^ — cj)), and vr are to be regarded as functions of {Jr, J^, 9r)] there is no 9^ dependence because 
of axisymmetry. We use the fact that vr = dR/dt = {dR/d9R)QR to obtain 
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We obtain in a similar way the equations 


_d 

dt 


d 


_d 

dt 


- y>) 

dvR 


d^4, 

dJ,, 


2J0 dR 


dJ,, 


— _ 


^ i?3 

3J2 


1 

^R 


^ch 


J ^ 
'W 


dVtE 


R4 


V^'{R) 


dR 1 

r p 1 

dJu ^R 

- VoiR) 


d^R 

dJ,, ■ 


(A25) 


(A26) 


(A27) 


Here is 1 for j/ = ^ and 0 for 12 = R. The set of three equations 1A25II through can be added to the other s 

to be integrated for an orbit, and they provide the additional values needed to evaluate the partial derivatives 1A24II . 
Initial values are dvRjdJ^ = d{9^ — (AijdJv = 0 at = t = 0 where R = i?min because vr = 9^ — cj) = Q there for all 
orbits. However the initial i?min values change with the actions, and initial values for the derivatives of R with respect 
to the actions are 

r D 1 dS 0_ U { u‘2 o T \ 

(A28) 
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They are obtained by differentiating the zeroth order energy equation. 

LAGRANGIAN DESCRIPTION, AND BOUNDARY INTEGRALS AS FLUXES 

The matrix analysis of ^uses an Eulerian description of phase space. Kalnajs (1977) gives an alternative Lagrangian 
description. An advantage of using a Lagrangian description of phase space is that it automatically includes any 
contributions from the motion of boundaries in phase space. Kalnajs’s Lagrangian analysis, with our definition m 
of Fourier coefficients, yields the formula 


Mjk{m,uj) = —47r' 
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with the unperturbed DF /o undifferentiated. Our equation ED follows after integrating m by parts with respect 
to Jr and J^. Those integrations introduce boundary terms at any boundary of the region of integration in action 
space unless one or other of /o and the term in square brackets vanish at that boundary. In the case of the wholly 
prograde DF of equation it^ . integration by parts with respect to gives precisely the boundary integral M® as 
defined in m, plus the area integral (|13 after also integrating with respect to Jr. This shows that the boundary 
integrals, which arise from the step function in the Eulerian description of a can be explained as due to the motion 
of that boundary. The boundary integral M® arises from the perturbation of a non-zero population of radial orbits. 

Lynden-Bell & Kalnajs (1972) had earlier used a Lagrangian analysis to calculate the perturbed angular momentum 
£2 and derived the result 
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They then integrated by parts to get an expression for the contribution L{hi,h2) to £2 from stars with angular 
momenta in the range {hi, /i2). Their result is 
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(B3) 


This formula of an area integral and two boundary integrals is the same as one obtains from applying the Eulerian 
equation IIA6II to the DF /o(Tfi, J^)H{J^ — hi)H{h2 — J^) which represents the stars with angular momenta in the range 
{hi, /i2). Lynden-Bell & Kalnajs interpret the boundary integrals as representing fluxes through the two boundaries. 
That interpretation, together with equations (PI) and m which give a physical significance to the real and imaginary 
parts of the matrix M, shows that neglecting the boundary integral terms i?7ll for a prograde disk of stars means 
neglecting the contributions to the total potential energy and angular momentum which arise from the perturbation 
of radial orbits. The numerical results reported in shows that this neglect can cause large errors. 










